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ABSTRACT: 



The two- and three-dimensional unsteady Navier-Stokes equations are 
solved numerically for the flow field about an impulsively started flat 
plate. In attempting to obtain an exact time dependent solution, several 
significant results were observed. First, with regard to the formulation 
of the differential equations themselves, it appears that Poisson's 
equation for the pressure field is a fundamental equation in as much as it 
allows us to solve for pressure most exactly at any given time. Secondly, 
the difference equations must be carefully and consistently formulated. 

In this research, a non-uniform lateral grid, a unique interpretation of 
the continuity equation, and "leap-frog" integration in time proved to be 
valuable techniques in obtaining an exact solution. 
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NOMENCLATURE 



At = Finite difference time step 

Ax = Finite difference grid spacing in x 

Ay = Finite difference grid spacing in y 

Az = Finite difference grid spacing in z 

T) = Transformed lateral coordinate = e ^ 

f = Arbitrary dummy dependent variable 

L = Characteristic length 

p, = Coefficient of molecular viscosity 

v = Kinematic viscosity 

p = Pressure 

p = Density 

Re = Reynolds number = pVL/p, 

t = T ime 

u = Velocity in x-direction 

v = Velocity in y-direction 

V = Characteristic velocity 

w = Velocity in z-direction 

x = Axial coordinate 

X = Body force vector 

y = Lateral coordinate 

z = Cross -streamwise coordinate 

u> = Successive over-relaxation parameter 

Subscripts 

i,j = Indicates components of a vector 
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NUMERICAL EXPERIMENTS WITH THE TWO -AND THREE-DIMENSIONAL 



UNSTEADY NAVIER STOKES EQUATIONS 
Gustave Hokenson 

INTRODUCTION 

Since Crocco (Ref. l) suggested seeking asymptotic solutions to 
approximate forms of the parabolic unsteady Navier-Stokes equations, in 
lieu of attacking the elliptic steady equations, many researchers (Refs. 2, 
3 , and 4 ) have successfully applied the technique to solve problems ranging 
from real gas nozzle flow to the hypersonic blunt body problem. The 
attempts generally are aimed at retrieving only the steady solution and 
this limited goal allows for extensive approximation within the differ- 
ential and difference equations as long as the steady flow is exactly 
represented by the equations. This approximation of the unsteady equations 
facilitates the numerical solution while supposedly not interfering with 
the uniqueness of the asymptotic solution. This technique, however, 
obviously negates the intermediate solutions from being actual representa- 
tions of the time dependent flow field. 

The essence of the numerical technique is to specify an "arbitrary" 
initial flow field, calculate the appropriate spacial derivatives in the 
Navier-Stokes equations, equate their sum to the time derivative of the 
respective dependent variable and step forward in time. Results of this 
investigation revealed that the combination of "arbitrary" initial condi- 
tions with the equation approximation previously mentioned creates some 
numerical problems which must be cloaked in various numerical subterfuges 
such as filters and arbitrary intermediate specifications of the velocity 
field. This is not meant to slight the technique for it does provide an 
easy method of obtaining apparently unique solutions for the steady 
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flowfield. However it was discovered in this research that in order to 
guarantee an accurate time dependent solution , instantaneously consistent 
with the equations and boundary conditions, much more care is needed both 
in the formulation and the numerical implementation of the differential 
and difference equations. 

The purpose of this investigation was to obtain the exact time 
dependent solution for the flow about a finite, infinitely thin flat plate 
impulsively started in its own plane into a uniform incompressible flow. 
Both the two- and three-dimensional Navier -Stokes equations were applied 
to the problem and the difficulties characteristic of both approaches are 
presented. A variety of techniques, all consistent with the equations and 
the boundary conditions, were employed to obtain exact time dependent solu- 
tions without recourse to non-rigorous numerical manipulations. 



DIFFERENTIAL EQUATION FORMULATION 



The Navier -Stokes equation for the flow of an incompressible 



Newtonian fluid can be written in Cartesian tensor notation: 




(i) 




( 2 ) 



where 
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j ax. 



(3) 
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With the aid of Equation 1, Equation 2 may be written in conservative form 
which is most amenable to exact numerical formulation (Ref. 5). In conser- 
vative form, Equation 2 becomes 



du. 

(u.u.) = X. - =• 
St Sx. v i j' i 



dP 



p dx. 



_2 

v V u. 



(4) 



Equations 4 and 1 form a closed system for the solution of the 
dependent variables u_^ and p since, with the appropriate boundary and 
initial conditions, there are as many equations as there are unknowns. 

These equations have been applied directly (Ref. 6) to obtain steady solu- 
tions to the Navier-Stokes equations with one of the momentum equations 
serving as an equation for pressure. Early in the course of the investi- 
gation a variety of problems were discovered with this approach simply in 
retrieving steady solutions and its utilization to obtain the desired time 
dependent solutions proved completely unfeasible. In lieu of solving a 
momentum equation for pressure, Equation 4 may be differentiated and summed 
to obtain a time independent Poisson's equation for pressure in terms of 
the instantaneous flow field. If we differentiate Equation 4 with respect 
to x^ and apply the Einstein slimming convention we get 




+ _L_ r JL_ 

ax Lax. 

J 





i a 2 p , 5 r 

p ax j ax i ax i L v 



2 1 
V u.J 



(5) 



Assuming the functions to be continuous, we can interchange the order of 
differentiation in the temporal convective and diffusive terms and obtain 
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(6) 
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for gravitational body forces we can write 




( 7 ) 



and for incompressible flow, Equation 1 states 




(8) 



Applying these conditions to Equation 6, we obtain the desired Poisson's 
equation for pressure 



o *2 

A 2 P ? 5 u i u i 

- - = V 2 P = - p 



dx.dx. 

l l 



dx.Bx . 



(9) 



With Equations 1, 4, and 9 we have a system of equations greater in number 
than the unknowns and we must choose the most fundamental set to solve. 

The specification of the apparent optimum choice was one of the goals of 
this research and is discussed in the next section. 

For generality in problem solving, we non-dimensionalize Equations 1, 

4, and 9 and for simplicity allow the new variables to have the same symbol 
as the terms in the original equations. It is understood from this point 
on that only the non-dimensional quantities will be discussed. The variable 
transformation is defined by 



u . 

l 



u . 

1 



V 






With this transformation the governing equations can be written 



( 11 ) 



du. 

1 

Bx. 



- 0 



3u 



1,3/ \ „ BP , 1 _2 

tt - + r (u.U.) = X. - + — V U. 

3t 3x. i j i Bx. Re i 
J i 



( 12 ) 



.2 

0 d u.u. 

V 2 P ^ 

Bx.Bx. 

i <3 



(13) 



where 



Re ^ 
M- 



For future use these equations are expanded for the three-dimensional case 
in rectangular Cartesian coordinates as follows 



3u | 3v [ 3w _ 
Bx By dz 



(14) 



du , a / 2 N , a / .x , b / x bp , i ra 2 u , a 2 u , a 2 ui /nc x 

Bt + Bx u + By (uv) + ai (uw) = - 3^ + Ri-Lrr + -2 + —2j (15) 

dx 3y oz 



|v + i. (uT)+ | (v 2 )+ A (TO) = .|z + E i[a! i + ^ + §] (16) 



dy 



dz 



dy Re 



'dx By Bz 



2 2 2 

dw , 3 / x d / x d / 2x BP , 1 Td w 3 w 3 w"| n _x 

at + SJ + a? < w > + ai S + 5? + ^2 + (17) 



o p p p p p p p p 2 2 2 

b p + il + Li - r a u . a v a w ~i ra_uv a vw a uw ] . 8) 

2 + 2 2 - ~ L 2 2 2j‘" LBxBy Bydz BxBzJ ^ ; 

dx By Bz Bx By Bz 



5 



For the two-dimensional case we can simplify these equations by dropping all 
terms involving w and z and their derivatives. In this case, the equations 
become 



du + dv 
dx dy 



- 0 



(19) 



du | d 
dt dx 



(u 2 ) + 




dP + _i_ |"dju + a 2 u 'j 

dx Re L 2 2.J 

dx dy 



( 20 ) 



dv 

dt 



+ 




+ 





_1_ !-afv 
Re L 3x 2 



+ 




( 21 ) 



J 2 . -.2 2 2 2 ,2 2 

d P + d P _ _ d u + p d uv + o_v_ ] 

-x 2 2 L -n 2 dxdy ^ 2 j 

dx dy dx J dy 



( 22 ) 



Just as the conservative form of the convective derivative is that 
which is most appropriate for numerical integration, the diffusive .terms 
also can be reformulated to most exactly satisfy the equations when finite 
differenced (Ref. 5). The reformulation of the diffusive terms is accom- 
plished by the following substitution in the three-dimensional equations 






2 2 
d v _ d w 

dxdy ” dxdz 



2 2 
d u d w 

dxdy dydz 



2 2 
d u _ d v 

dxdz dydz 



(23) 



(24) 



( 25 ) 
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and for the two-dimensional system 



a 2 u _ a 2 v 
ax 2 " ‘ 9x3y 



(26) 



a 2 v 



3y 



5 2 u 

SxSy 



(27) 



where the appropriate continuity equation has been differentiated to obtain 
Equations 23-27. After performing these substitutions the 3-d equations 
become 



Su [ 5v | Sw _ 

ax ay az 



( 28 ) 



3u , 3 / 2\ B \ \ dP , 1 P3 P- 

at + to (u > * to (u,,) + to (uu) ' ‘ to + ite Lrr 



2 2 2 2 

ap ^ i ra u ^ a u _ a v _ a w 

2 axay axaz_ 



Sy 5 z 



(29) 
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rr + — (uw) + — 

at ax ay 



2 2 2 

a p . a p , a p 

,2 2 2 

ax ay az 



( i 3 / \ SP , 1 , 

(v ) + _ (w) , . _ + _ 


ra 2 v 


a 2 u 




axay 


(vw) + (w 2 ) = - 

az ' ' az Re j 




.2 
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. 2 ' , 
Sy 


_-2 2 ,2 2 2 2 n 
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2 

ra UV 


>2 

! a vw 


■ L 2 \ 2 V 2j 2 

ax ay az 


_axay 


ayaz 



2 2 

a w , a v 



ayaz 3z 2 



: J (30) 



a 2 u 

axaz ayaz 



a 2 v^ 



j (3D 



(32) 



For the two-dimensional case all w and z terms and their derivatives may be 
dropped and we retain the following equations. 



a_u 

ax 




(33) 
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( 34 ) 




(35) 




(36) 



DIFFERENTIAL EQUATION APPLICATION 



As was mentioned earlier, there is a choice with regard to the 
differential equations which can be used to solve the problem. As a general 
rule, it was found necessary for both the continuity equation and Poisson T s 
equation for pressure to be satisfied at each time step, including t = 0 , 
to guarantee an accurate time dependent solution. Based on this principle 
each streamwise velocity component was calculated from its respective 
momentum equation, the cross -streamwise velocity from the continuing equa- 
tion and pressure from Poisson’s equation. If there is no distinct cross - 
streamwise direction, for example in a wake or cavity recirculation region, 
a momentum equation is solved for each velocity component, Poisson’s equa- 
tion is solved for the pressure, and continuity is checked and guaranteed by 
the use of correction factors (Ref. 7) in Poisson’s equation in the terms 



which cannot be arbitrarily set equal to zero if continuity is not being 
used to solve for one of the dependent variables. 

For the particular problem we have attacked the differential equation 
application is as follows 
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u from u momentum 



2- d v from continuity (37) 

P from Poisson f s 

u from u momentum 

w from w momentum 

3- d (38) 

v from continuity 

P from Poisson's 

In order to obtain exact time dependent solutions, it was found 
necessary to specify physically meaningful and compatible initial condi- 
tions throughout the flow field, not "arbitrary initial conditions" as is 
so often quoted in the literature. The compatible and physically meaning- 
ful initial conditions for this investigation were specified as follows 

u = 0 on the surface of the plate 

u = 1.0 everywhere except on the plate 

2- d (39) 

v = calculated from continuity 

P = calculated from Poisson's 

u = 0 on the plate 

u = 1.0 everywhere except on the plate 

3- d v = calculated from 2-d continuity (40) 

w = calculated from 3-d continuity 

P = calculated from Poisson's 

There may be some question about the approach for the 3-d problem, however 
no inconsistent results were found using this approach. It appears that the 
calculation of w from 3-d continuity after calculating v from 2-d continuity 
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would result in a null matrix for w. That this is not the case is observ- 



able from the way in which the differential continuity equations were 
implemented, using locally spacial averaged values of u, v, and w. This 
procedure is outlined in the next section and the appropriate difference 
equations are listed in the Appendices. 

FINITE DIFFERENCE EQUATION FORMULATION 
Based on the approach outlined in the previous section, the spacial 
partial derivatives of the dependent variables in each of the differential 
equations were finite differenced using three point central difference 
formulas. Special care was taken to avoid known inaccurate difference 
models (e.g. the Du-Fort-Frankel representation of diffusive terms) when we 
are seeking exact time dependent solutions. Substitution of the appropri- 
ate difference expressions into the momentum equations allows for the 
calculation of the instantaneous temporal derivatives of the dependent 
variables u and w based on knowledge of the current dependent variables . 

The differencing procedure for the momentum equations is outlined in 
Appendix A, followed by a listing of the respective calculations of the 
temporal derivatives as they appear in the computer program. In both the 
two and three-dimensional situations Poisson f s equation was used for the 
calculation of pressure. Numerical experiments were conducted using the 
y-momentum equation for the solution of P, however, since we dealt with the 
exact equations the time derivative of v must be formed and added to the 
appropriate spacial derivatives in order to calculate the partial deriva- 
tive of P with respect to y. With the pressure on the freestream boundary 
known, it was in theory possible to integrate the equation inward from the 
outer boundary. However, numerical instabilities arose from the inability 
to calculate physically compatible initial conditions since at this time v 
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is known only at one time. In addition, this method of calculating P has 
inherent accuracy limitations and this study revealed the necessity of 
calculating P to a very high degree of accuracy, including the initial 
conditions . This discovery served as the impetus for using Poisson's equa- 
tion which was solved using classical successive overrelaxation techniques 
with all derivatives in the equation being represented by three point 
central difference models. The finite differencing of Poisson's equation 
is presented in detail in Appendices B and C along with a listing of the 
encoded form as it appears in the full computer program. 

The most crucial finite differencing formulation which was encountered 
was that of representing the continuity equation by a formula which is both 
accurate and physically meaningful. The method which was settled upon is 
equivalent to a third order Runge-Kutta approach and the details of its 
formulation along with the computer listing is presented in Appendices B 
and C. 

In each of the equations which were finite differenced the cross - 
streamwise partial derivatives were formulated with the option of using a 
non-uniform grid spacing. This allows for the most efficient use of 
computer memory and concentrates the calculation grid in the regions where 
the dependent variables vary most rapidly. The central difference formulas 
for both the first and second order partial derivatives with a non-uniform 
grid are developed and presented in Appendix A. 

FINITE DIFFERENCE EQUATION APPLICATION 

The finite difference equations discussed in the previous section were 
applied to the solution of the incompressible flow about an impulsively 
started flat plate in the manner outlined in Equations 37 > 38? 39? and 
The physically compatible initial condition field was found and the 
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appropriate equations were used to continue the velocity field downstream 
in time with Poisson's equation being solved for the new pressure field. 

When explicitly updating the velocity field in time, a stable solution is 
obtainable for a limited range of time step increments. Karplus (Ref. 8) 
obtained numerically stable solutions with a similar system of equations 
with a restriction for the time step which can be written in non-dimensional 
variables 



A 7 A Z Re ^2 ^2 



( 41 ) 



Equation 4l proved to be a completely satisfactory stability criterion in 
this work. 

In an effort to obtain the most exact solution (to second order in At) 
the mechanism of using effective central differencing in time was introduced 
and is included as an option in the computer program. In general, an 
independent variable f can be calculated at a future time from a second 
order Taylor series expansion in time as follows 



f(x,y,z,t + At) 



f(x,y,z,t) + |£ (x,y ,z,t) (At) 




(x,y,z,t) 



iAtf 

2 



Since f(x,y,z,t) is known and — can be calculated from the differential 

ot p 

of 

equations, it only remains to calculate — -. For the Navier-Stokes equa- 

af at d 

tions — is an explicit function of the dependent variables and their 
ot o 

d^f 

spacial derivatives and it therefore follows that — p may be calculated by 

at* 

differentiating the governing equations with respect to time and inter- 
changing spacial and temporal derivatives. However the implementation of 
these equations is lengthy and the need to solve a second Poisson's equation 
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S P 

(for rr) increased the computing time intolerably. Because of these 
ob 

factors , it was decided to attempt to obtain the same second order accuracy 
using the basic difference equations appropriate to a first order approach. 
This was accomplished by calculating the dependent variables at time t + At/2 
with forward finite differencing and using these values to calculate the 
spacial derivatives in the Navier-Stokes equations at that time. The 
dependent variables axe then stepped forward from time t to time t + At 
based on the temporal derivatives evaluated at t + At/2. The implementa- 
tion of this approach is most clearly understood from the flow chart of the 
computer program presented in Appendices D, E. The inherent susceptibility to 
numerical instabilities which this "leap-frog" technique exhibits should be 
damped by the strongly dissipative terms in the equations (Refs. 9> 10 , and 
11 ). 



Finally, with respect to obtaining the solution throughout the flow 
field, some technique must be specified to meet the boundary conditions. 

Apart from setting u=v=w=0on the surface of the plate, some procedure 
is necessary for setting u, v, w, and p on the boundaries of the calculation 
region. Two methods are available and have been proven to be essentially 
equivalent. First, the values of the dependent variables which have been 
calculated in the interior of the region may be extrapolated to the 
boundaries using the arguments of symmetry and uniformity across the 
boundaries to generate the appropriate extrapolating function. Second, the 
equations themselves may be applied to the boundaries if special considera- 
tions are taken which give expression to those values needed outside of the 
calculation region, based again on symmetry and uniformity arguments. The 
latter method was chosen as the preferable from the standpoint of consistency 
of accuracy and error generation throughout the flow field. 
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COMPUTER PROGRAM 



The computer programs which encoded the solution of the finite 
difference equations in Appendices B and C are listed in Appendices D and E, 
along with the flow charts appropriate to each. With regard to programming 
the system of finite difference equations and several mundane considerations 
such as the memory capacity of the computer and the accuracy required of the 
system vs. what the equations are demanding of it are deserving of some 
mention. First, with regard to the 3-d. program (with the inherent 3-d 
arrays of dependent variables), is the question of computer core storage 
capability. The minimum number of 3-d arrays which must be dimensioned is 
equal to the number of dependent variables. Because the exact difference 
equations are being used, the calculations of new variables must be made in 
conjunction with the storage of old ones and it is apparently not possible 
to utilize only the minimum number of arrays. However in this program the 
temporal derivatives are written on discs immediately after being calculated, 
each as a dummy scalar. When the calculation of the dependent variables at 
a future time (either t + At or t + At/2) is being initiated, the appropriate 
discs are rewound and the scalar strings are read in sequence to provide the 
temporal derivative of the proper dependent variable at that point in space. 

In addition to the artifice of using virtual memory in the form of 
simple disc writing, an additional tack was tried. This was to map the 
infinite space in y into a finite space with the transformation 

7) = 

if this is done, the appropriate partial derivatives become 
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By 





-if 4 

BTl 




and the lateral distances can be stretched such that there are fewer points 
laterally but there is denser data in the region of steepest gradients . 

This transformation also allows the external boundary conditions to be set 
somewhat more exactly and, in conjunction with the non-uniform grid size as 
established for y, should give the optimum use of total lateral space per 
point of computation. Because the results were adequate at the Reynolds 
numbers which we were studying, the transformed coordinate was considered an 
unnecessary complication for this problem. It is felt however that many of 
the computational difficulties associated with the exact numerical solution 
of lower Reynolds number problems may be alleviated by the use of this 
artifice . 

With regard to the questions about the accuracy of the solutions, two 
specific areas of investigation were dealt with in this research. First was 
the utility of computing in double precision (real*8) and second was the use 
of higher than second order approximations to the spacial partial derivatives. 
The computations were all performed with a variety of grid sizes and accurate 
solutions were obtained with spacial and temporal grid sizes large enough to 
allow significant differences to be calculated in single precision (real*4) . 

If the non-uniform lateral grid is set up with extremely small spacing near 
the plate, the ensuing small differences may force the use of double 
precision to retain the most accurate solution. In general it is felt that 
this is the only situation which should force one to use double precision. 
Experimentation with the use of higher order approximations to the 
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spacial partial derivatives was carried out with the following two representa- 
tions of the spacial derivatives 



df 



~ f .i+2 



+ 8 Vi - 8 

12 Ay 



’ + f 

' 3-1 . 1-2 




- f .i+2 + 16 f .i+l - 3° fl-1 - 

12 fiy 2 



The results using these models were identical to those using the second order 
approximations with adequate grid size and the additional complication of 
increasing the thickness of the boundary on which the boundary conditions are 
to be set could be avoided by using the standard formulas. Therefore, neither 
the use of double precision nor the inplementation of "more accurate" finite 
difference models was found to be of general use. 



RESUITS AND DISCUSSION 

The exact unsteady two and three-dimensional Navier-Stokes equations 
have been applied to the solution of the incompressible flow about an 
impulsively started flat plate. Some typical results are presented in 
Figure 2 which gives the transient development of both the downstream and up- 
stream flow on the centerline and in the plane of the plate. For our purpose, 
it is not particularly startling that we can solve this specific problem but, 
more importantly, that we can numerically solve for the time dependent flow 
field using the two and three-dimensional Navier-Stokes equations. In so 
doing a set of non-rigorous rules has been obtained which can be used as 
guidelines in the application of the equations toward solution of more 
complex problems . 

The Navier-Stokes equations written in strict conservative form and 
physical variables can be integrated to obtain accurate time dependent 
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solutions providing care is taken in the formulation of the problem. At 
all times during the calculation, the continuity equation and Poisson's 
equation for pressure must be satisfied to ensure a valid solution while 
each streamwise velocity component must satisfy the appropriate momentum 
equation. In addition, the initial conditions must satisfy the continuity 
equation and Poisson's equation for pressure in order to allow a "natural" 
flow development. The boundary conditions can be satisfied either through 
application of the equations to the boundary with the use of symmetry and 
uniformity arguments or through the use of polynomial extrapolation from the 
interior calculation region. 

In essence, Poisson's equation for pressure has been introduced as a 
fundamental equation in the numerical solution of the full equations 
because of the need for generating a compatible initial pressure field with 
a time independent equation and also because of the need for calculating the 
pressure field to an extremely accurate degree. This accurate pressure field 
calculation and the specification of physically compatible initial conditions 
eliminated the need for all numerical artifices such as filters to cover such 
purely numerical phenomena as velocity overshoot . 

Second in importance to the adequate solution of Poisson's equation is 
the satisfaction of continuity equation at each point in time and space. 
Apparently, only one finite difference formulation of this equation, that 
presented in the Appendices, is entirely adequate in obtaining exact time 
dependent solutions of the problem. 

The governing equations were finite differenced using a central 
difference approximation for all spacial derivatives. The DuFort-Frankel 
formulation of the diffusive terms was avoided in order to guarantee the 
most exact representation of the spacial derivatives at each instant in time. 
In order to obtain apparent central differencing in time the dependent 
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variables at time t + At were calculated using the spacial derivatives 
calculated at t + At/2 where the t + A/2 values were obtained by pure 
forward differencing. To obtain stable magnitudes for these time steps , the 
approximate equation given by Karplus was used as a guide once accurate 
spacial grid sizes were chosen. 

Beyond the formulation of the numerical problem, the operational 
problem of obtaining sufficient computer memory space was the most difficult 
to overcome. For both the two and three-dimensional problems, the use of 
disc storage was applied to help alleviate the problem. As is clear from 
the flow diagram of the program, the intermediate calculations were stored 
on discs as dummy scalars and re-read as matrices when the dependent vari- 
ables were being updated. This technique allows the use of a minimum number 
of dimensional arrays. In addition, the use of a non-uniform grid and 
mapping the calculation region into a finite space were applied to help 
alleviate the space problem and optimize the number of calculation points at 
those places where they are needed. 

Several further investigations are now underway to determine the effect 
of the initial conditions and the equation approximation on the uniqueness 
of the asympototic solution. In addition the two and three-dimensional 
instability problem is being studied from the standpoint of comparing the 
effect of two and three-dimensional disturbances and finite vs. infinitesimal 
disturbances on the classical instability maps. Finally the full compressible 
problem is being formulated in a manner identical to that presented here. 
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CALCULATION REGION FOR FLOW OVER A FINITE FLAT PLATE. 




FIG. 2 

TIME DEPENDENT SOLUTION FOR THE FLOW OVER 
AN IMPULSIVELY STARTED FINITE FLAT PLATE. 
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APPENDIX A 



Finite Difference Formulas for Non-Uniform Spacing 




2 

We let F = Ay + By + C and expand around y . to find A , B , and C . 

0 

F = A Dely 2 (j) + B Dely (j) + C 

F. = C 
J 

F = A Dely 2 (j-l) + B Dely (j-1) + C 

J — 

Solving for A and B we get: 



Dely (j) + Dely (J-l) 

= Dely (J) Dely (J-l) 
Dely (J) + Dely (J-l) 



F - F 

f.i+1 j 


F . .. - F. ... 

+ J_1 , J 1 


L Dely ( J) 


Dely (J-l)J 


F - F 

j .1+1 3 


F. , - F. 
J-l J 


^Dely 2 (J) 


Dely 2 (J-l) 



Now that A, B, 



and C is known. 



we can differentiate F at 



y = y . and it is 



clear that 



SF 

By 



= B 




= 2A 
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APPENDIX B 



Finite Difference Formulation of 3-d Navier-Stokes Equations 



I. u Momentum Equation 



A 



B 



C D 



E 



F 



G 



H 



2 





A=(U( 1 + 1, J,K)*t 2-U( 1-1, J,K)**2) /( 2.0E0*DELX) 

B 1 =U ( I » J* 1 » K ) * V ( I » J +1 t K)/ ( D EL Y ( J )**2 ) 

B2 = U( I , J,K)*V(I , J t K)M 1.0E0/DELY( J-l ) F * 2-1 . OEO/DELY ( J ) M 2 ) 

B3=U ( I t J 1 " 1 1 K ) *V ( I , J-l» K )/DELY( J-l)*=»2 

B=DELY( J)*DELY( J-l )*( B1+B2-B3 )/ ( DELY ( J ) + DEL Y ( J -1 ) ) 

C=(U( I, J,K + 1)*W( I , J ,K + 1)-U( I , J, K-l)~ W(I , J, K-l) ) / ( 2.0E0* OELZ) 
D=(PHI ( I+1,J»K)— PHI ( I- 1 » J » K ) )/( 2.0 EO fe DE L X ) 

E2A=V( I +1 , J+l ,K) /DELY( J) 2 

E2B=V ( 1+1, J,K )*( 1.0E0/DELY( J-l)** 2- 1. OE 0/DE LY( J)**2) 

E 2C =V ( I + 1 » J-l » K ) /DELY (J-l ) **2 

E2=DELY( J)*DELY< J-1)F ( E 2A+ E2B-E2 C ) / ( DEL Y ( J ) +DELY ( J-l ) ) 

E 1 A=V ( 1-1, J+1,K)/DELY(J )**2 

E1B=V( 1-1 , J,K) K 1 .0E0/DELY( J-l ) **2-1.05 0/ DELY ( J ) **2 ) 

E 1C=V( 1-1, J-l, K) /DELY( J-l)** 2 

E1 = DELY ( J ) *DELY (J-l )* (E1A+E1B-E 1C ) / ( DEL Y ( J ) +DEL Y( J-l ) ) 
E=(E2-E1)/(2.0F0*DELX) 

F=(W( 1+1, J,K+1)-W( 1+1, J,K-1)-W( I-1,J,K+1)+W(I-1,J,K-1) ) / 

1 (4.0E0*DELX*DELZ) 

G1 = U( I , J+l ,K) /DE LY ( J) 

G2=U ( I , J, K) M 1.0EO/DELY( J-1H-1. OE 0/DEL Y( J ) ) 

3 3 = U( I , J-l t K)/DELY( J-l ) 

G=2.0E0*( G1-G2+G3) /( DEL Y( J ) +DE LY( J-l ) ) 

H=(U( I, J,K+1)-2.QE0*U(I, J,K)+U( I, J,K-1) ) /DELZ**2 
UTEE=- ( A+B+C+D+REI*(E+F-G-H) ) 
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II. w Momentum Equation 



A B C D 



E 



F G H 



2 

Sw r duw Bvw + dw 

Bt _ L Bx By Bz 



+ — + i ■{ ^ a + ^- [ — I 

Bz Re [BxBz Bz LByJ 



2 2 
B w B w\ 1 

. 2 ' . 2JJ 
Bx By 



A=( U( 1+1, J *K)»W( 1+1 t J»K)-U( 1-1, J ,K)*W( 1-1* J»K ) )/ ( 2.0E0*DELX) 
3 1 = V( If J+1,K)*W(I , J+1,K)/DELY(J)**2 

B2=V< I, J,K)*W( It J, K)*( 1.0E0/DELY( J- 1 ) * * 2- 1 . OE O/DE L Y ( J) **2) 

B 3=V ( I f J-l.K) *W( I , J-i f K)/DELY( J-l )*»2 

B=DELY( J-1)*DELY( J )*( B 1+B 2-B2J / ( DELY( J) + DELY ( J-l ) ) 

C = (W( I, J, K+l) **2-W( I , J ,K-1) **2) /< 2.0E0*DELZ) 

D= (PHI ( I, Jf K + D-PHI ( I, J,K-1) )/( 2 .OEO^DELZ ) 

E=(U( I + 1 .J t K+l ) — U ( I +1 »J*K-1 )-U( 1-1 , J»K+1 )+U (1-1, J,K-1) )/ 

1 ( 4.0E0*DELX*-DELZ ) 

F2 A=V ( I * J+l » K +1 J/DELY ( J )**2 

F2B = V( I »J»K+1)M1.0E0/DELY(J-1)** 2-1 .OEO/DELY(J )**2 ) 

F2C=V ( I,J-1,K+1)/DELY(J-1)**2 

F2=DELY(J)* DELY (J-l)* ( F2A+F2B-F2C )/ ( DEL Y ( J )+DELY( J-l) ) 

F 1A=V( I » J+l » K-l ) /DE LY( J) **2 

F 1B=V ( I t J , K-l ) * ( 1.0/DEL Y( J-l )*$ 2- 1.0/DE LY( J)**2) 

F 1C = V ( I * J— 1 » K — 1 ) /DELY ( J-l ) **2 

Fl=DELY(J) sfc DELY( J - 1 ) M F 1 A F 1B-F 1C ) /< DE L Y ( J ) +DELY ( J-l ) ) 
F=(F2-F1)/(2.0E0*DELZ) 

G = ( W( I + lfJtK)— W( 1-1 ♦ J ♦ K) ) /<2.0E0*DELX) 

H1=W( I » J + 1 » K ) /DEL Y( J ) 

H2=W( I, J, K )*( 1.0 EO/ DELY (J-1J+1.0E0/ DEL Y( J ) ) 

H3= W( I , J-1,K)/DELY( J-l) 

H=2.0EO*(H1-H2 + H3)/(DELY( J J+DEL Y( J-l) ) 

WTEE=-( A+B+C+D+REI* (E+F-G-H) ) 
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III. Poisson's Equation for Pressure 



2 2 2 

BP | o P BP 

2 2 2 

Bx^ By 3z^ 



B 


C 


ra 2 u 2 




Ux 2 


ay 2 


sfp 


+ <L| + 


Bx 2 


ay 



,2 2 

+ + 
3z 2 





E 


F 


G 


/a_ 


-auv I 


+ A. n 


_ 2 

'Bvw , ( d uwl | 


lax 


L By J 


3z L 


. ByJ BxBzjJ 



be written 



2 2 2 

^ = A1 + A2 + A3 - P(l,J,K)*H 

Bx^ By Bz^ 



combining equations and we can solve for P 

P(I,J,K) = (A + B + C + D + 2.0 (E + F + G))/H 



where A = A1 + A2 + A3 and H are defined in the subsequent listing. Since 
equation is not an explicit solution for the pressure field, we must employ 
some iterative process such as successive over relaxation which is used here. 
With this approach, equation is written 

P(I,J,K) = (1 - a>) P( I , J,K) +oo (A + B + C+ D + 2.0 (E+F + G))/H 



A1=(PHI ( I+1»J»K)+PHI ( I-l»JfK)) / DELX* *2 
A2=(PHI (I»J»K + 1)+PHI( I, J.K-in /DELZW 

A3=2 . 0 EO* (PHI (I , J + 1,K)/DELY( J»+PHI( I, J-l, K)/0ELY( J-l) )/ 
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1 ( DEL V ( J ) +DEL Y( J-l) ) 

A= A1+A2 +A3 

B = (U(I + ltJ»K)< l < e 2-2.0E0* t U(I , J » K ) j?* 2+U ( I - 1 , J , K ) ** 2 ) / DELX**2 
C1=V( I t J+1»K)**2/DELY( J ) 

C2=(V( I ,J,K)**2)*(1.0E0/DELY( J- 1 ) +1.0/DELYC J) ) 

C 3=V( I » J-l ,K) **2/DELY( J-l) 

O2.0E0MC1-C2+C3) /(DELY( J )+DEL Y( J-l) ) 

D=(W(I » J,K+1)**2-2.0E0$W( I» J,K)**2+W( I, J ♦ K- 1 )#*2 ) /DE|_Z*v2 
E2A=U( I + lt J+l» K) *V(I + 1, J+l »K) /DELY< J)**2 

E2B=U( I +1, J,K )*V ( 1+1, J,K)*( 1.0E0/DELY( J-l )**2-l. OE O/DE L Y ( J ) ** 2) 

E2C=U( 1+1 , J-l ,K) *V( 1+1 , J-l ,K)/DELY (J-l ) ** -2 

E2= DEL Y ( J ) *DEL Y( J-l)* ( E 2A + E 2B-E 2C ) / < DELY ( J ) +DELY ( J-l ) ) 

El A=U ( I-l,J+l,K)*V(I-l,J+l,K)/DELY(J)*y2 

E1B*U( I -It J,K) > V( I-l» J,K) * ( l.OEO/DELY (J-l )**»2-l . 0 EO/ DEL Y (J )*% 2 ) 

E1C=U( I-1,J-1,K)*V< I- l t J-ltK) /DELY( J-l) **2 

E1 = DELY( J)*DELY( J-l )* ( E 1A+ E 1B-E 1 C ) / ( DEL Y ( J )+DELY( J-l) ) 

E = ( E 2-E 1 ) /( 2-OEO^DELX) 

F2A=V( I , J+1,K+1)*W( I, J + l t K+1)/DELY< J)# » 2 

F 2B = V ( I , J,K+1 ) W C I , J,K+1)*<1 .OEO/DELY( J-l ) *'*2-1 . 0 E 0/ DEL Y ( J ) ^^ 2) 

F2C=V( I »J-1»K+1)*W( I »J-1»K+1) /DE LY( J-l ) ** 2 

F2 = DELY ( J )*DELY ( J-i )M F2A+F 2B-F 2C ) / ( DEL Y ( J ) +DEL Y( J- 1 ) ) 

F 1 A=V( I » J+ 1 » K— 1 ) *W(I»J+1,K-1)/DELY(J)**2 

F1B=V( I,J,K-1)*W(I,J,K-1)M l.OEO/DELY (J-l )**2-l . OEO/ DEL Y ( J ) * * 2 ) 
F1C=V( I ,J-1,K-1)*W(I,J-1,K-1)/DELY(J-1)**2 
F1 = DELY( J) *-DELY( J- 1) F ( F 1 A+F 1 B-F 1 C ) / ( DEL Y ( J ) +DELY ( J-l ) ) 
F=(F2-F1)/(2.0E0$DELZ) 

G2=U(I«-1, J,K+1)*W( 1+1, J,K + 1 )+U( I-1,J»K-1 )*W(I-1t J»K-1) 

G 1=U( 1-1, JtK+1) *W( I -It J,K+1)+U( 1 + 1, J,K-1 )*W(I+1 , J »K-1) 
G=(G2-G1)/(4.0E0*DELX*DELZ) 

H = (2. 0E0/DELX**2+2.0E0/DELX**2+2.0E0/ (DELY( J)*DELY(J-1) ) ) 

APHI=( 1. OEO- OMEGA )vPHI ( I , J , K ) + 0 MEGA*- ( A+ B + C+ D+2 . 0 EO* (E+F+G) )/H 
RESID=A6S(APHI-PHI(I, J,K) ) 

ERR=MA X 1 ( ERR »RESI D) 

PHKItJ ,K)=APHI 
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IV. 3-d Continuity Equation to Solve for w 



dw fdu + 5v ~j 
5z " Ldx dyj 



UEX2 =U ( I + ltJ + l*K)+U{I + ltJ-l»K)+U( I + ltJ + l»K + l)+U( I + 1»J-1,K+1) 

U EX 1 = U { I-l»J + l t K)+U{ I~1 * J — 1 . K ) +U C I — 1 » J + lfK+1) +U( I— 1»J- 1»K+1) 
UEX = ( UEX2-UEX1)/{8.0E0*DELX) ’ ' 

V 1= ( V { 1+1 ,J+1,K) +VC 1+1, J+1,K+1) +V(I-1, J+ 1,K)+V(I-1 , J+l , K+l ) )/ 



V2=( VI ! 1+1, JtK) *V( I + 1 , J ,K+1) +V( 1-1 , J,K) + V ( 1-1, J,K + 1 ) )/4.0E0 
V3={V(I+i»J-l»K)+V{I+ltJ-l»K+l)+V{I-ltJ-ltK) + V(I-lJj^lt K+I ) 
1 4 . 0 EO 

VY1=V1/DELY( J) t*2 

V Y 2 = V 2 * { 1.0E0/D5LY{ J)**2-l. OEO/DELYI J-l )** 2) 

VY3 = V3 / DELY ( J— 1 )**2 

VWY=DtL Y( J) 1 DELY{J— l)4 t (VYl+VY2-VY3)/(OELY(J)+DELY(J-l) ) 

W( I.J,K+l)=W(I f J,K)-DELZ*(UEX+VWY) 



)/ 



To evaluate ^ at I,J,K + l/2 the following method was employed. First v 
o z 

a t y , y., y. was found by averaging the four v values in the x-z plane 
J -L J J 

at a given y. With these average values, the non-uniform grid derivative 

Sv / 

formula was applied to form — at I,J,K + 1/2. In a similar manner, an 

average u at I + 1 and I - 1 was found by averaging the four values in the 

y-z plane at each I. These average u T s were then central differenced to 

form r— at I,J,K + l/2. Finally, using these derivatives evaluated at 
dx 

I,J,K + l/2, the w field is calculated using Equation using a central 
difference method. 
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V. 2-d Continuity Used to Solve for v 



dv _ du 

dy " dx 

The identical procedure was applied to the finite difference solution of 
this equation in two dimensions as was outlined for the 3-d continuity 
equation. The principle difference being that the appropriate averaging 
of the dependent variables now occurs on lines instead of planes in the 
3-d matrix of grid points. 

C1=U( 1+1, J+lrK)+U(I + l*J*K) 

C2=U(I-1, J+1»K)+U( 1-1, J,K) 

V( I , J+l ,K) =V( I , J,K>-DELY(J)*(C1-C2)/(4.0EO+DELX) 

VI. 3-d Continuity Used to Solve for v. Following the procedure out- 
lined in Section IV of this Appendix, we can interpret the three-dimensional 
continuity equation to be an equation for v, here again using averaged 
values of the dependent variables in the difference formulation. 



U EX2= ( U ( 1 + 1, J,K-1)+U( 1 + 1, J+ 1,K-1 ) +U( 1 + 1, J,K+1)+U( 1 + 1, J+l.K+l) )/ 
14. OE 0 

UEX1= < U(I-1, J,K-1) +U( 1-1, J+1,K-1) + U(I-1 , J,K+1)+U( 1-1 ,J + 1 ,K + 1 ) )/ 
14.0E0 

UE X=( UEX2-UEX1 )/(2.0E0^DELX) 

WZ2=(W( I+l,J t K+l) +W( 1+1, J+1,K+1)+W(I-1, J>K+1)+W( 1-1 » J+l » K+l ) )/ 
14.0E0 

WZ1 = ( W( I + lt Jf K-1)+W(I+1, J+lfK-l)+W(I-l, J,K-1)+W( 1-1, J + 1,K-1) )/ 
14.0E0 

WZE={ WZ2-WZ1 J/(2 .0E0*DELZ) 

V( I » J+ 1 » K ) =V( I , J,K)-DELY( J)^(UEX + WZE) 
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APPENDIX C 



Finite Differencing of 2-d Navier-Stokes Equations 



I. u Momentum Equation 



A 



B C 



D 



E 




A=(U( I+L, J )**2-U( I-l» J )**2)/(2.0E0*DELX) 

81 = U( I » J+l ) *V( I i J + l )/DELY(J )**2 

B2=U( I, J )*V( I, J ) M 1.0E0/DELY( J-l) **2-1. OEO/DELY( J)**2) 
B3 = U( I, J-1)*V(I, J-l )/DELY( J-l)+*2 

B =OEL Y( J ) *DE L Y( J - 1 ) M B1+B2-B3) / ( DELY( J ) +DELY { J-l ) ) 

C = ( PH I ( I + ltJ)-PHI( I — 1 •» J I )/(2.0E0*DELX) 

D2A=V ( I +1 , J+L )/DELY ( J )**2 

D2B = Y{ I+1,J)*( 1 . OEO/DEL Y( J-l)t*2-1.0E0/DELY(J)* fe 2) 
Q2C=V ( I+1.J-1)/DELY(J— 1)**2 

D2=DELY(J)*0ELY(J-1 )* ( D2A+D2B-D2C )/ ( DELY ( J )+DELY( J-l ) ) 
D 1A=V ( 1-1, J+l )/DELY( J )**2 

D1B=V( 1-1, J)+(l .0E0/DELY( J-1)**2-1.0E0/DELY(J )**2) 

D 1C = V( I-l»J-l)/DELY{J-l) k ? t 2 

D1 = DEL Y ( J )*DELY( J-l )t ( D1A + D 1B-D 1C ) / { DEL Y ( J ) +DEL Y ( J-l ) ) 
D = (D2-D1)/(2.0E0*DELX) 

E 1 = U( I ♦ J+ 1 ) /DELY ( J ) 

E2=U ( I , J )*( 1.0E0/DELY( J-l ) + 1 . OE 0 /DEL Y( J ) ) 

E 3=U( I , J-1)/DELY( J-l ) 

E=2. OE 0* ( E1-E2+E3) /{DELY( J)+DELY( J-l) ) 

UTEE=-( A+B+C+REI*(D-E ) ) 



II. Poisson's Equation for Pressure 



B 



C 



D 




> p 

From finite differencing 2—— + - — — , we can show that 

By 
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= A1 + A2 - P(l,J)*H 

3x By- 



Substituting this into equation, we obtain 



P(I,J) = (A + B + C + 2.0 D)/H 



To obtain a form appropriate to successive over -relaxation, we can write 



P(I,J) = (1 - to) P(I,J) + co(A + B + C + 2.0 * D)/H 



A1=(PHI (I+1,J)+PHI( I-l.J) )/DELX**2 

A2=2.0EO*(PHI (I, J + 1)/DELY(J)+PHI( I,J-1)/DELY( J-l) )/ 

1(DELY( J)+DELY(J-1) ) 

A= A1+A2 

B=(U( 1+1, J)**2-2 .0E0*U( I, J )**2+U< 1-1, J )**2) /DEIX**2 
C1=V( I , J+1)**2/DELY( J) 

C2=V ( I , J)**2*( 1.0E0/0ELY( J- 1 )+ 1 . 0E0/DEL Y( J) ) 

C3 = V( I , J-1)**2/DELY (J-l J 
C=2.0E0*(C1-C2+C3) / ( DEL Y( J ) +DE L Y ( J-l ) ) 

D2A=U ( I +1, J+l )*V< 1+1, J + 1)/DELY( J)**2 

D2B = U(I +1 » J ) * V( 1+1 , J)*(1.0E0/DELY ( J-l )**2-l .OEO/DELY( J )**2) 
D2C=U ( 1+1, J-l^Vt 1+1, J-1)/DELY( J-l)**2 

D2=DELY ( J)*DELY( J-l)* ( D2A+D2B-D2C ) / ( DEL Y( J )+0ELY( J-l)) 

D 1A =U( I - 1 , J+1)*V(I-1, J+l) /DELY( J)**2 

D1B=U ( I — 1 ».J ) *V( I-ltJ)^(1.0E0/DELY(J-l)) ti 2-1.0E0/DELY(J)) t ^2) 
D1C=U ( I -1 t J-l ) * V ( 1-1, J-1)/DELY( J-l)**2 

D1=DELY( J )*DELY( J-1)M D1A+D1B-D1C ) /( DELY ( J)+DELY( J-l) ) 
D=(D2-D1)/ (2.0EO*DELX) 

H=2.0E0/DELX**2+2.0E0/(DELY(J)*DELY( J-l ) ) 

APHI = ( 1 . 0E0-3ME GA )*PHI ( I , J ) +0ME GA* ( A+B + C + 2. 0E0* D ) /H 
RESID=ABS( APHI-PHI ( I, J ) ) 

ERR=MAX1( ERR f RE SI D) 

PHI ( I , J ) =APHI 



III. 2-d Continuity, Solving for v 



3v 3u 
By - 3x 



To form — at I.J + l/2, average values of u at I + 1 
By 

by averaging the values of u along appropriate J, J + 
differencing of these average values yields at I, 



and I - 1 were found 
1 lines . Finite 
J + 1/2. 



UEX2=( U(I+lt J)+U( I+lt J+l) )/2.0E0 
UEX1=( U(I-1, J)+U( 1-1, J+l) )/2.0E0 
UEX=( UEX2-UEXl)/( 2.0E0*DELX) 

V (I t J+l ) = V( If J)-DELY( J )*UEX 
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APPENDIX D 

I. FLOW CHART FOR 3-D PROBLEM 
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APPEND I X D 



II. LISTING 0= COMPUTER PROGRAM FOR 3-0 CASE 

THE FOLLOWING FORTRAN CODED FINITE DIFFERENCE COMPUTER PROGRAM 
SOLVES THE FULL TIME DEPENDENT NAVIER STOKES EQUATIONS FOR THE 
THREE DIMENSIONAL FLOW ABOUT A FINITE, INFINITELY THIN FLAT PLATE 
IMPULSIVELY STARTED IN ITS OwN PLANE. "LEAP-FROG TIME-WISE INTEG- 
RATION IS INCLUDED AS AN OPTION. EXTENSIVE DISC WRITING IS EMPLOY- 
ED TO SAVE COMPUTER CORE SPACE. ADDITIONALLY, THE EQUATIONS 
THEMSELVES ARE APPLIED TO T H F BOUNDARIES OF THE CALCULATION REGION 
IN LIEU OF EXTRAPOLATING THE VARIABLES CALCULATED WITHIN THE 
REGION TO THE BOUNDARIES. 



THE FOLLOWING PARAMETERS MUST BE SPECIFIED 
ICAIC LE ME NE LP1 LP2 NP1 MP2 PE DEL X OELY DELZ OMEGA ERRTOL D E 



Q 



v 



> ,■« >■ 4 f 4 s:. n f ■. 4 ** i- f '■< se V. , "• •k % ft* :* * f y e. k. /. . »; * v. .* * ,f 

f- 7 

* ICALC SPECIFES THE USE OF "LEAP-FROG" IN T * 

t. 'i 

* LE IS THE LENGTH OF THE REGION IN DEL X STEPS 

r> 

* ME IS the HEIGHT OF THE REGION IN DELY STEPS * 
£ NE IS THE WIDTH OF THE REGION IN DELZ STEPS *: 

..c 

i L PI IS THE START OF THE PLATE 

#/ 

* LP2 IS THE END OF THE PLATE 
$ 

NP1 IS ONE SIDE OF THE PLATE 
V NP2 IS (JNE SIDE OF THE PLATE 

* RE IS THE REYNOLDS NUMBER BASED ON U AND L 

i 

a DELX IS THE GRID SPACING IN X 

t DELY IS THE GRID SPACING IN Y 

* DELZ IS THE GRID SPACING IN Z 

> OMEGA IS THE RELAXATION PARAMETER IN PGI SSONS 
•;* ERRTCL IS THF ERRCR TOLERANCE IN PCI SSONS 

* PE IS THE FR EEST REAM PRESSURE IN PSF 
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c 

c 

r 

c 



c 

c 

c 

c 



r 



c 



t 

c 

r 

C 



C 

r 



QE IS THE FREESTREAM DYNAMIC PRFSSURE IN PSF 
4 NOTE THAT NE MUST 6F AN ODD INTEGER 

: NCTE THAT DELY 

ALL VAR I ABLE S ARE 

7 

1 * *»• «*- '• . *. i .L i ’J r « Jf / * */ L 4 . 

0 [ M ENS I CN U ( + 0 , 30 , 2 i ) , V ( AO , 

D I '•“> E S I ON DELYI30) 

L f i =L E-l 
L r 2=LF - 2 
L F3 = LE- 3 

M 2 1=NE- 1 
MF 2=^ c -2 
ME3 =ME-3 

N E1=NE- 1 
NE2=NE-2 
N E3=NE- 3 

NPM=(NE+1 ) /2 
NPtfl=NPM- 1 
N d Y2 = N p M + 1 

RFI=1.0c0/R p 

PHI E=PE /( 2. 050 QE) 

A = 1 .0 EC/DELX 
B=I.OEO/DELY( 1) 

C=1.0CO/rELZ 

DELT = i.0E0/(A+0.1E0' B+O.l B0*C+2 .0 EOtRR I * ( A ,: : >2+B') * 2+C : '' 2) ) 

C?LTl = Ccl T 
DELT2=D2LT1/2.0E0 

IC = 0 

INITIALIZE THE CALCULATICN FIELD 
DO 2 I = 1 , L E 
DC 2 J = 1 , ME 
DC 2 K = 1 » N E 



■i 

IS an ARRAY DEL Y (ME ) 

DIMENSIONED FILE, ME, NF) * 

V 

I v Jf \ A * *l 3i - •* r , .a.**’ «- ,*n 

30, 21), W(40, 30, 21) ,«>HI ( 40,30,21) 



3 ^ 



U( It J,K ) = 1 .OEO 

V ( I » J » K ) =0 .0 EO 
W( I ,J,K)=O.OEO 
PHI ( I , J * K ) = PHI E 

2 CONTINUE 
C 

C SET PLATE TO ZERO VELOCITY 

DO 3 I = LP1 *.LP2 
DO 3 K=NP1,NP2 
U ( 1,1 , K ) = 0 .OEO 

3 CONTINUE 
C 

C CALCULATE LNITIAL V WITHIN THE REGION FROM 2-D CONTINUITY EO 

DO 4 J = 1 » ME 1 
DO 4 1 = 2, LEI 
DO 4 K=1 , NE 

C1=U( I+1,J+1,K)+U(I+1, J,K) 

C2=U(I-1»J+1,K)+U(I-1»J,K) 

V( I , J+l ,K) =V( I , J , KJ-DELY ( J )- (C1-C2 )/ ( 4 .0 EO^DELX ) 
m CONTINUE 
C 

C SET UPSTREAM V l-IELD 

DO 5 J = 1, M E 
DO 5 K= 1 , N E 
V( 1, J ,K ) =0. OE 0 

5 CONTINUE 
C 

C EXTRAPOLATE DOWNSTREAM V FIELD 

DC 6 J = 1 , M E 
DO 6 K=1 , NE 

VUE, J ,K)=3.0E0? V( LEI, J , K ) - 3. OE 0* V< LE2, J,K)+V(LE3, J,K) 

6 CONTINUE 
C 

C CALCULATE W FIELD WITHIN CALCULATION REGION FROM 3-D CONTINUITY EO 

CG 7 1=2, LEI 
DO 7 J =2 , ME 1 
DO 7 K= NP M, NE 1 

U£X2=U( 1 + 1 , J+1,K) +U( 1+1 ,J-1 ,K)+U( 1 + 1, J+1,K+1)+U( 1 + 1, J- I,K+1) 
UEX1=U( I-1,J+1,K ) + U( I- I , J- 1 ,K) + U( 1-1 , J+l , K+l) +U( I -1 , J-l ,K+1 ) 

U FX= ( UEX2— UEX1 ) / ( 8 . OEO--‘DEL X ) 

VI = < V< I +1 , J+l ,K) +V ( 1+1 , J+l * K+l ) +V( 1-1 f J + l, K )+V ( 1-1, J + l, K + l ) ) / 

14. OEO 

V2 = (V( 1+1, J,K )+V( I +1 , J,K+1 )+V( I - 1, J,K )+V( I- 1, J,K + 1) ) /4. OF 0 
V3 = < V( 1+1 , J-l ,K) +V( 1 + 1 , J-l , K+l) + V( 1-1 , J-l ,K)+V ( 1-1 , J-l, K+l ) )/ 

14. OEO 

VY1 = V1/DELY( JH*2 

VY2=V2M 1. OEO/DELY( J) **2-1. OEO/DELY( J-l )*' 2 ) 

VY3=V 3/ DE LY ( J- 1 ) : +2 
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V wY = CEL Y ( J ) DE LY{ J-l) M VY 1+ VY2- VY3 ) / { DE L Y ( J ) + DE L Y ( J-l) ) 
w(I,J,K + l> =WU, J,K)-DELZ* (UEX+VwY) 

7 CONTINUE 

C 

C CALCULATE W IN THE PLANE OF THE PLATE 

DC 8 1 = 2, LEI 
DC 3 K= NPM, NE1 
J = 1 

UEX2=U( 1 + 1, J + l, K ) +U( 1 + 1, J+ I ,K) +U( 1+ 1, J<- 1 , K+ 1)+U( I +1 , J+l , K+l ) 
UEX1=U( 1-1 ,J + 1 ,K)+U( I-1,J + 1,K)+U( I-1,J + 1,K + 1)+U( I- 1 , J+1,K+1) 

U C X = ( UE X2-UEX1) / ( 8.0E0* DELX) 

V 1= (V ( I + 1»J+1»K.)+V ( I + 1,J+1,K+1) + V( I— 1»J+1,K) + V( I— 1 , J+l ,K+1) ) / 
14.0E0 

V2=(V( I+l,J f K)+V(I+I,J ,K+1)+V{I-1,J,K)+V( I-ltJfK+1) ) / 4 . 0 EO 
V3 = ( V ( I+1,J + 1,K )+V ( I + 1,J + 1,K+1) + V ( I-1,J + 1,K) + V( I — 1, J+l, K+l) ) / 
14. DEC 

V 2=-V 3 

VY 1 =V 1 / DEL Y ( J )* +2 

V Y2 = V2 - ( 1. OFO/DELY{ J) W-l.OtO/ DELY( J) • ' 2 ) 

VY3=V3/CELY{ J K : 2 

VhY = OELY{J ) D F L Y { J ) . ( VY 1+VY 2-V Y 3 ) / < DEL Y { J ) +DEL Y( J ) ) 
w( I,J,K+1)=W(I , J,K)- r ELZ* (UEX+VWY) 

8 CONTINUE 
C 

C CALCULATE W AT THE TOP OF THE CALCULATION REGION 

DO 9 1 = 2, LEI 
DC 9 K=NPM,NE1 
J = ME 

U EX2=U ( I +1 , J , K ) +U ( 1 + 1, J-1,K)'+U( I + 1,J»K+1)+U(I + 1,J-1,K+1) 
U r Xl=U( 1-1 , J ,K) + U( 1-1 , J-l , K) +U< I-1,J,K+1)+U(I-1,J-1,K+1) 

U E X- ( UE X2- U p XI ) / ( 8 . OEO- DELX) 

VI = (V{ I +1 , J,K)+V< I+i, J,K+1 )+V( 1-1, J,K )+V( I- 1, J,K + 1 ) ) / 

14. OEsO 

V2=(V( 1+1, J,K)+V< 1 + 1, J,K+1)-+V( 1-1, J,K)+V< I-1,J,K+1) ) /4.0E0 
V3 = (V( I +1 , J-l ,K) + V( 1+1 , J-1,K+1 ) +V( I- 1 , J-1,K )+V( I - 1, J-1,K+1) ) / 

1 4 . 0 E C 

V Y 1 = V 1/ CELY ( J >•* - 2 

V Y 2 = V 2 (l.GE0/DE.LY(J)*-'2-<-l .OEO/DELY(J-1 :*■ 2) 

VY3=V3/DELY( J-l)’ • 2 

VWY = DELY( J )-• DEL Y (J-l)* ( VY 1 + VY2-V Y 3 ) / ( DEL Y { J ) + DELY( J-l) ) 

W( I , J, K+l) =W( I , J,K) -DE LZ^ (UFX+VWY) 

9 COf-IT INUE 
C 

C INSERTION OF Lc c T HALF REGION W VALUES BY SYMMETRY 

DO 10 1=2, LEI 
DC 10 J =2 , ME 1 
DC 10 K =NPM2» NE 
K 1 = N E + 1 — K 
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*< I,J,K1) = W! I, J,K) 

10 CCNTINUE 
C 

C CALCULATE U ON THE FRONT OF THE CALCULATION REGION 

DC 11 J=!, ME 
DO 11 K = 1 » NE 
W(1,J,K)=0.0E0 

11 CCNTINUE 
C 

C CALCULATE W AT THE EXIT 3 C THE CALCULATION REGION 

D n 12 J =1 » ME 
CO 12 K=1,NE 

w(LE,J,K)=3.0E0. : W!LEl,J,K)-3.0E0 r ‘W{L E2, J,K)+W(LE3,J»K) 

12 CONTINUE 
C 

C SET W ON THE PLATE SURFACE EQUAL TO ZFRO 

DO 13 I=LP1,LP2 
CO 13 K=NP1,NP2 
W(I ,1,K)=0.0E0 

13 CONTINUE 
C 

C CALCULATE INITIAL P FIELD 

r 

WRITE(6,810) 

£10 FORMAT! 1H0,20X,' I NTO IC POISSONS') 

GO TO 200 
150 IC=1 

WRITE! 6,811) 

611 FORMAT ( iH0,20X, 'OUT CF IC POISSONS') 

I T ER=1 
1 A CONTINUE 

IF( ICALC.EQ.O) &0 TO 18 
I TE R=0 
REWIND 16 
REWIND 18 
C 

C. STORE PREVIOUS FIELD ON TAPES 

W R I TE ( 1 6 ) U 
WRITE! 18) W 
C 

18 CONTINUE 
C 

REWIND 19 
REWIND 20 

C CALCULATION WITHIN THE REGION 



CO 


19 


1=2, LEI 


DO 


19 


J =2 , ME 1 


DO 


19 


K = 2 , ME 1 
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C :.'-LC<JL£- T IC\l CF 1JT F-ry 313 MQy: fv 7 uM - QliAT I C fv 

-- (J ( i +1, J , K )- •• 2-U ( 1-1, J, K i » 2)/(2.0E0"DELX) 
ji=U( 1 » J+l »K)*’ V( I , J + l ,K) / ( DELY ( J )-•' 2 ) 

S 2 =U( It J,KJ V( It J»KJ' ( 1.0E0/CELY( J-l) k 2 - 1 . OE 0/D E L Y ( J ) ' 2) 

E3 = U( I » J-l ,KJ* V ( 1 , J-l ,K )/DELY ( J- 1) 2 

5=0 EL Y( J) -OELY( J-l) (B1+B2-B3)/ ( CELY ( J ) +DELY ( J-l ) ) 

C=(G( 1 ,J,K + 1 ) W( I , J ,< + l)-U( I , J , K-l) . W( I , J ,K-1) ) / ( 2. O t O DtLZ) 
r =( PHI ( 1 + 1 , J,K)-PHI ( 1-1 , J,K ) )/ ( 2 .0 50‘DELX ) 

= 2/- = V( 1 + 1» J+ 1» K) /DEL Y( J) !■ 2 

- 2 *5=V ( I + 1 1 J t K ) ’ ( l.OEO/DEL Y( J— 1 > * 2-l.OF 0/DcLY( J) ' ’ 2) 

- 2 C=V( i + 1 , J-l ,K) /CELY( J-l )' " 2 

r z= DEL Y ( J ) L 1 L Y( J- 1 ) (E2A + E23-E2C)/(DELY(J) +DELY ( J-l J ) 

C 1A=V( I — 1 1 J + l ,K )/DELY( J )» 2 

E 1 B = V(I- 1 ,J,K) ( 1 . OrO/DELY ( J-l ) * • 2-1 . 0 EO/ C ELY ( J ) « -2 ) 

EiC=V( 1-1. J-1,*)/P5L Y( J-l) 2 

- 1=DFLY ( Ji • OE t Y { j-i ) ( clA+ c 18-t 1C )fl DEL Y( J ) +CELY ( J- 1 ) ) 

E - ( C 2- '£ 1 ) / 1 2 . 0:. 0 Dr LX) 

P=U( 1+1, J, *+!)-*( i+ 1 , J,K-l)-*( 1-1, J,K+1) + U(I-1, J,K-1) ) / 

1 ( 8 . 0 r 0 Of LX DtLZ) 

G 1= U( I , J+l.K) /DEL Y ( J) 

GZ = U ( I » J » K ) ( 1 .0 =0/ CELY (J-l H-l.Ot O/DEL Y( J ) ) 

G ?•- U( I , J-l » K ) /Z\ LY ( J-l ) 

C--2.QEG' (G1-G2+G3}/CELY( JJ+DELY( J-l) ) 
h= cue 1 * J»K +1 J-2.0E0 J( I,J,K)+'J( I, J,K- 1 ) 1/DELZ^ 2 
U T 5r=-( A+B+C+D+Rcl (E+F-G-H) ) 
w- IT E ( 19,903) UT EE 
VU- Ft FMAT( 1 PE 14. 7) 

r. 

C CALCULATION OF WT FPC M 30 MOMENTUM EQUATION 

A=( U( 1 + 1 , J » *< ) 5 w( 1+1 , J, KhUlI-I » J .K) >W ( I -It J,K ) ) / ( 2.0EC- DEL X) 
3 1= V ( I » J + 1 » K ) * r( I »J-*-l»K)/QZLY( J ) 1 v 2 

62 = V( I , J,K ) W ( I , J, K) v ( 1. CEO/ DEL Y (J-l)’ i 2- 1 . CE 0/D E L Y( J ) ^ ' 2 ) 
o3 = V(I tJ-1 t K ) » W( I ,J-l,K)/DELY(J-l)>-2 
B = T EL Y ( J- 1 ) DEL Y ( J ) ( B 1+6 2- B 3) / ( DEL Y( J ) +DE L Y ( J-l ) ) 

C = (W(I , J,K + l)^2-h( I,J,K-1 )*-’2)/(2.0E0*DELZ) 

D = (PHI ( I , J t< + l) - p H! ( J , J»K-1 ) ) /( 2 .OcO*DSLZ) 

r =(U( I + ltJ»K + l)-U( I+1»J»K— 1)-U( I— 1»J»K+1)+U(I— 1»J»K— 1) 1 / 

1 ( 4. OEO • DE LX ; DELZ ) 

F2A = V( I , J + i.K + l ) /PrLY( J)> 2 

- 2 E = V( 1 , J, K+l ) ( 1 .OcO/DELY( J-l)’ 2- 1 . OE 0/ D EL Y ( J ) * 2 ) 

C ^C=V( I ,J-i » K+ 1 ) /CELY( J-l) ‘ *2 

F 2 =DELY( J ) • P = LY( J-l)* ( F2A + F 26-F2C) / (DcLY( J)+DELY( J-l ) ) 

F 1 A =V ( I , J+1,K-1)/DELY(J)» ‘2 

F 1 0 =V ( I , J , K- 1 ) ; ( 1. C/DELY( J-l ) ‘ * 2-1 .0/0 ELY ( J ) **'2 ) 

F1C=V ( 1 , J-1,K-I)/DELY( J-l) *-’2 

C x-DELY(J) DELY(J-l)* (P1A + F16-F1C)/ (DELY (J )+DELY( J-l ) ) 
F=(F2-Fl)/(2.0t0 DtLZ) 

G= (W ( I +1 , J, K ) -w ( 1-1, J , K ) )/ ( 2-.0E0- DELX) 
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9 04 

19 

C 

c 



c 

c 



c 

c 



H 1 - W ( i , J+1,K ) / DEL y { J ) 

H2=W(J , J,K). (1.0 50/DFLY(J-1 ) +1 . 0 E0/ DELY ( J ) ) 

H2 = W( I , J-1,K) /0ELY( J-l ) 

H = 2 .0 EG- ( H1-H2 + H3 )/ ( DELY ( J ) +DEL Y( J-l) ) 

WTE E =-( A + B+C + D+PEH (E + F-G-H) ) 

W R ITE ( 20, 904 ) WTEE 
FORMAT ( 1PE14.7 ) 

CONTINUE 

CALCULATION ON THE SIDES 
DO 20 I =2, LEI 
DC 20 J = 2 » M E 1 
K = 1 

CALCULATION OF UT FROM 3D MOMENTUM EQUATION 

K-l TERMS HAVE BEEN ADJUSTED BY UNIFORMITY ARGUMENT 

A = ( U( 1 + 1, J *K) *»2-U( 1-1 * J,K) 2) / (2.0E0 - QELX ) 

B1.*U( I , J + l. K )* V( I » J + 1,K )/( DELY( J)**2) 

32 = U(I , J,KK« V (I , J* K) M L.OEO/DELY (J-l )• ' 2-1 .OEO/DELYI J)' -2) 

B 3 = t)( I * J-l, KM V( I * J-l ,K) /DFLY( J-l )~ 2 
B = CELY ( J )» DELY( J-l)4 (B1+B2-B3) /( DELY( J)+DELY( J-l) ) 

C = (U(I* J,K+1M W(I, J,K+1)-U( I,J,KH-W(I,J,K) ) / ( 2r0 EO* DEL Z ) 

D = ( PH I i 1 + 1*J*K)-PHI(I-1*J,K) )/( 2.0E0/DELX) 

E2A=V( I +1 » J + l » K )/DELY(J)*f'2 

E 2B = V( I +1 , J,K) M 1. OEO /DELY (J-I 2-1 .OEO/C ELY (J')^"2 ) 

E2C=V ( 1 + 1. J-1,K)/DELY( J - 1 ) ^ - 2 

E2 = DELY ( J) ^DELY ( J-l )*■ ( E2A + E2 6-E2 C ) / ( DEL Y ( J ) +DEL Y( J- 1) ) 

E1A=V( 1-1, J+1,K) /DEL Y( J) * *- 2 

E 1 B=V ( I - 1 * J * K ) ; ( 1.0E0/DELY(J-i) 2- 1 . OE O/DE LY ( J) * ► 2) 

E 1 C- = V ( I -1 , J-l ,K) / DELY (J-l )*«2 

E 1= DEL Y ( J ) -DEL Y( J-IL- (E 1A + E 1B-E 1C) / (DELY ( J) +DELY (J-l) ) 

E= ( E2- El )/ (2 .0E04DELX ) 

F = ( W( 1 + 1* J,K+1) -W(I+1 , J,K) - W ( I - 1 , J*K+1 )+W ( 1-1 , J,K ) )/ 

1 (4.0E0 : DEL # DEL Z ) 

G1=U( I * J + l * K ) /DELY (J) 

G 2 = U( I , J,lOM 1. OE O/DE LY( J-l ) + l . OEO/CELY ( J ) ) 

G 3 = U ( I , J-1,K )/DELY(J-l) 

G = 2a OE 0* (GL.-G2 + G3 ) / ( DFLY( J ) +€'ELY ( J-l ) ) 

H = <U( I , J ,K+l)-2. OE O' U( T * J * K ) +U( I , J * K) ) / DE LZ^ 2 
UTEE=-( A+B+C + D + REI-^( E + F-G-H) ) 

WRITE (19* 903) UTEE 

CALCULATION OF WT FROM 3D MOMENTUM EQUAT ION 

A = ( U( I + 1 * J * K ) *W( I + 1*J,K)-U(I-1 , J , K ) ’ W ( I -1 *J,K) )/ (2 .0 EG- DE LX ) 
B1*V( I , J + 1,K )*W< I , J + 1,K)/DELY( J )«■« 2 

B2*V(I , J,K) *W( I . J, K)M1 .OEO/DELY ( J- 1 ) * f 2- 1 . OEO/ DEL Y( J ) * "' 2) 
B2=V( I * J-1,KM W( I , J-1,K) /DELY( J-l)% > 2 
B = CELY(J-1 ) * D E L Y ( J )" ( B 1 +B2- B3 ) / ( DEL Y( J ) +DEL Y( J- 1 ) ) 

C=(W(I , J,K+1)*< 2-W( I , J » K) •“•'•2) / (2.0E0 1 DELZ) 
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D = ( PHI ( I , J ,K+1 ) -PHI ( I , J, K ) )/ (2 .0 E0 ' DELZ ) 
E=(U(I+1,J,K+1)-U(I + 1,J,K)-U(I-1,J,K+1>.+ U<I-1,J,K))/ 

1 (4 .OEO DELX* DELZ ) 

F 2 A = V ( I ,J + 1,K + 1)/CELY(JP -2 

F2E=V( I , J,K + 1P ( l.OEO/DELYt J- 1 ) * - 2- 1 . OE 0/ CE LY ( J ) * * 2 ) 

F2C = V ( I , J-l , K + l )/DFLY ( J-l ) < J 2 

F 2 = 05 L Y ( J ) jELY(J-I) ( F2 A+ F 2 d-F 2 C ) / ( DEL Y ( J ) +DELY ( J-l ) ) 
F1A=V( I , J + 1,K )/D=LY( J i -‘ v 2 

F 1 b = V ( I , J , K) ’• (1 . 0/ DELY ( J-l )* 2-1 .O/DELY ( J 2 ) 

F1C = V( I , J-ltK) /DE LY( J-l) ‘ ’2 

F1 = DELY < J ) 'GELYl J-l ) * ( F1A + F 1B-F 1C )/( CELY( J)+DFLY( J-l) ) 

F = (F2-F1) /(2. OEO *• DELZ) 

G= (W( 1 + 1, J,K )-w( I- 1, J,K) ) /( 2.0E0- DELX) 

HI =w ( I, J+l.K ) /DEL Y ( J ) 

H2 = w( I , J , K) ‘ ( 1 . OE O/DE LY ( J-l ) +1 . 0 FO/ DELY ( J ) ) 
ri .;■= w ( I, J-l, K) /DEL v( J-l) 

r\-tL .OEO (H1-H2 + HJ) / ( D E L Y ( J ) +DEL Y( J-l ) ) 

V. i EE --( A+fi+C+D+Pr I (E+F— G— H) ) 

*■- ITT (20,904) 

2 0 C r M: MJ E 

r 

l C «LCU L AT I C • ) UN THE: SIDES 

C Ki-1 7 c r K S HAVE BEEN ADJUSTED BY UNIFORMITY ARGUMENT 

0 21 1=2, LEI 

Cl 21 J = 2 , M 2 1 

K = 

C CALCULATION QF JT FROM 3D ^OMENTUM EQUATION 

A=(U(I+1,J,K)» ’ 2— U( I — 1 , J , K ) - 2 ) / ( 2 .GEO" DEL X ) 

61 = U( I , J+1,K)«V( I , J+1,K) /(DELY( J)*» 2) 

b2=U(I»J»KHV(I»J,K) ( l.OEO/CcL Y( J-l)'- - 2-1. OEO/DELY( J) • -2) 

5 3 = *J( I , J-l , K) (■ V ( I , J-l , K ) / DELY ( J-l )v"2 

6= DEL > ( J ) FDEl Y( J- 1 ) ( B 1+B2-B3) /( DELY( J) +DELY( J-l ) ) 
C=<U(I,J,K)W(I,J,K)-U(I,J,K-1) -W(I,J,<-1))/(2.0E0-CELZ) 

C = < PHI ( 1+1 , J ,K) -PHI ( I— 1 , J,K) )/ (2. OEO- DELX ) 

F 2A = V ( 1+1, J + 1,K)/0ELY( J)f >2 

E? r : = V ( 1+1 , J,K) . ( 1 ,O c O/DELY ( J-l ) ' '2-1 .OFO/DELY( J )“ 2) 

:-2f=V( I + ltJ— 1,K) / D E L Y ( J-l) ! ,l 2 

F2=DFL.Y ( J ) '-DELY ( J — 1 ) ' ( E2A+E 2B-E ?C )/( DEI Y( J)+DELY( J-i) ) 

E i A »V ( I — 1 , J+ 1 , K ) / D E L Y ( J ) o - 2 

E1B = V( I - 1 , JtK) ■ ( 1. CE O/CE L Y( J-l) ’ 2-1 . OE O/DE LY ( J ) * - 2 1 
El C=V ( I— 1 »J — 1,K )/DELY( J— 1 ) *= 2 

El = Df LY (J) ' DEL Y( J-l )*■ (E1A+F1E-E1C)/ (DELY ( J ) + DEL Y ( J-l ) ) 
E=(E2-El)/(2.0t0'DFLX) 

F=(W( 1+1, J,K)-W( 1+1, J,K-n-W(I-l, J,K)+rt( I-1,J,K-1) ) r 
1(4. OEO DELX' DELZ) 

G! = U< I , J + 1 , K ) /DEL Y ( J ) 

G2=U( I , J, K) . ( 1.0 =0/DELY (J-l 1+1 .OEO/DbLYU ) ) 

6 -) = U( I , J- 1 , K ) /DELY( J-l) 
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G = 0 < G1-G2 + 3 3) / ( Z B L Y ( J ) + 0 E L Y ( J-l ) ) 

h= (u( It J, K )-2.0E 0 - U( I , J,« ) + tf( I , J ,K-1) ) /0ELZ- 2 

UT= = --( ^+E-»-C + U+K£l - ( t + F-G-H L) 

rt 1 ' I u f { 1 G » 903 ) UT . F 

C.-LCULATICK' CF WT FF'M 30 MOMENTUM C QJATICM 

A= ( J( I + lt J»* ) a( I + 1»J,K)-U(I-1»J,K) • W(I-itJ,iO)/(2.0E0 DELX) 
-1=V(I,J+1,K) '*(I,J+1,K)/0ELY(J) - 2 

•^2=V(I,J,K) V!(I»J,K) ( 1.0E0/dELY ( J- l ) 2 -1 . 0 =0 / 0 EL Y ( J ) * ■ 2 ) 

2i=V( I » J- 1 , K ) --W ( I , J-l, K) /DBLY( J- 1) 2 

3=r-.LY(J-l) DEL Y ( J ) { 6 1+62-S2 ) / ( CFL Y ( J ) +DEL Y( J- 1 ) ) 

C - ( i t J»K) • 2- w ( I , J , K- 1 ) " 2 ) / ( 2 . 0 FC ' 3 L Z ) 

. = (PHI ( I, J,K)-PHI ( It J»K-1) ) /( 2. OtO OcLZ ) 

;=<U(I + 1,J,K)-U(I-H,J,K-1)-U(I-1,J,K)+J(I-1,J,K-1))/ 

1( R.CJ^O'CELX CFLZ ) 

F2A = V( I , J-H,K )/DELY(J ) 2 

r 2e = V( I , J » K) { 1 . 0 F 0 / 1 >b LY ( J — 1 ) 2~1 .0 -0 / OtLY ( J ) 2) 

T ,J-I,K i/DELYi J-l) 2 

: -U' LY { J) OtfLYC J-J. ) (F2A + f-2c-FiC )/( OELV(J )+DELY( J-l )) 
r IA = V( I , J+ i-t K— i ) /£.: u Y t J-) * • 2 

Flr=V(I,J, K.-1 ) LI . 0 / EL Y ( J - i ) 1.0/0LLY( J ) ; A 2~) 

FlOiM! » J-l » K-l ) / LFLY ( J-l )’ 2 

- i- CF L Y { J ) " 0 E L Y ( J - 1 ) ( F 1A +• f 1--F 1 C ) / ( DE L Y ( J) +DE LY < J-l ) ) 

F = { Fz-Fl )/ (2 .O c O'"CcLz ) 

G=(W(1+1»J,K)-W( 1-1 , J, K) )/ (2.0tC- CELXi 
H l = to ( I , J + 1 . K ) /D5LY< J) 

H2=h( I , Jr K)" ( 1 .3 FO/TELY ( J - 11+1 . 0 F 0/ DEL Y t J ) ) 

H 3 = w{ I , J-l ,K) /CELY( J-l) 

H- ^ .0 F 0 ( H1-H2+H3)/(DFLY{ JJ+OEL Y( J-l) ) 
to7EF=-( A+B+C+D+REI ' (F+F-G-H) ) 
wFIT c ( 20,904) toTEF 
21 O.TIM'F 

CALCULATION ON THE BC.TTCM OF THE CALCULATION RE3IQi\ 

J-l T t c MS HAVE BEEN ADJUSTED ON SYMMETRY ARGUMENTS 
D" 22 1=2, LEI 
D: 22 K = 2 » N E 1 

J--1 

CALCULATION GF UT FROM 3C MOMENTUM EQUATION 

A = ( U ( I + 1 , J , K ) -r 2-U ( I - 1 , J , K ) * 2) /(2.0E0< DELX) 

=U( I , J + l ,K) i V( I , J+l , K ) / ( DELY ( J )A 2 ) 

3 2= U( I , J , K ) V(I, J,K)' ( 1 . 0 c 0 / C r L Y ( J— 1 ) ! 2~1 . GEO/DELY ( J ) *• 2) 

EO = U( I , J + I,x )' V( I, J + l ,K )/DlLY( J )**• V2 

d = F EL Y ( J ) -^ OF L Y ( J ) (3i + B2-B3)/(DFLY(J)+DELY(J)} 

r » ( U ( I, J » K +1 ) ► to { I,J,K + 1)-J( I , J , K— 1 ) * W ( I,J,K-1) )/(2.GE0 i DELZ ) 

: =(PHI { 1 + 1 , J , K ) - PH I (I-lfJ»K) )/ (2.G EO* DELX ) 

-2A = V ( I+i, J + 1,K )/DEL Y( J ) ' - 2 
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c / L- = tf ( 1 +1 , J> K » • ( i . OEO/DE L Y( J ) • 2- 1 . OF 0/ D E L Y ( J ) — 2 ) 

E2( =V( i+1, J+-l,M/DfcLY<J)- 2 

t 2C =- r 2C 

r 2 = Q E L Y ( J 1 - 'CEL Y( J )~ ; ( E24 + E 2 d-E2C ) /(CELY( J ) +DELY< J) ) 
l i*-=V( I -1 , J+l , K)/ DELY ( J ) ‘ > 2 

1 Ir =V< I - 1 , J , K) ( 1. OEC/fTELY ( J) • 2 -1 ,0 EG / 0 E LY ( J ) > 2 ) 

r xC=V ( I-l» J + l, K )/DtLY{ J )* ’2 
E 1C =-E 1 C. 

c'l=DFLY(J) s -DELY( J) ■ (E 1A+E1B-E1C) / (DELY( J)+DELY(J) ) 

E = ( 52- = l )/ (2 .060 ■' DELX ) 

F = ( w( 1 + 1 , J ,K+1 )-W( 1+1 , J, K-l )-W ( 1-1, J, K+l ) +W ( 1-1, J tK-1 ) ) / 

1(4. 0=0 DELX DELZ) 

G 1 = U ( I , J+l * K )/D5LY { J ) 

G 2 = 'J( I , J,K) ( 1.0E0/0ELY(J) + 1.0EG/DELY(J) ) 

G3 = U< I, J+1,K)/D=LY( J) 

G =2 .0E0- : (G1-G2 + G3 )/ (PELY(J ) +DEL Y ( J ) ) 

4 = ( U( I , J, K+l) -2. 0E0 : U( I , J, K) + U( I , J ,K-1 ) ) /CELZ*--2 
JTE= = -( A+B+C+D+REI' (E+F-G-H) ) 

.. = IT= (19,903) UTEE 

CALCULATION CF WT fFGK 30 MOMENT UK EQUATION 

A-(U(Ifl»J»K)< to ( I + 1 ♦ J » K ) — U ( I— 1 » J » K i‘ w ( I — 1 » J > K ) ) / ( 2.0E0-- DELX) 
B 1 = V( I,J+l.K)iw( I , J+1,K)/DELY( J)> '2 

6?=V( I » J ♦ K L W { I » J , K )' (1.0E0/DELY(J M * 2- 1 . OF C/DEL Y{ J ) *■ 2) 

B 2 = V( I , J+l ,K) ’ to( I , J+l, K) / OELY( J )* ‘Z 

3 "■ = - F 3 

5=DFLY( J) >-DELY( J) 5 (B1 +62-33 ) / ( C EL Y fJ ) +OEL Y ( J ) ) 

C = ( W( I , J, K+l) » *2-W(I , J , K- 1 ) ** 2) /(2.0E0*DELZ) 

D = (PHI ( I, J ,K + 1)-PHI ( I, J,K-1 ) )/( 2 . OEO/DE L Z ) 
c = <U( 1+1 , J,K+1) -U< 1+1 , J,K-1 )-U( I-WJ,K+1 )+U <1-1, J,K-1) )/ 
1(4.060- DELX' DELZ ) 

F2A = V( I , J + l ,K+1 )/DELY( J)* 2 

F 2? = V( I i J * K+ 1 ) ( l.GEO/DELY( J) s >• 2-1.0E0/DELY(J ) ■• 2 ) 

F/.C=V( I , J + i, K+nVCFLYtsT)* 2 
C ?C =-F 2 C 

F2 = DELY(J) OELV(J) ( F 2A+ F2 3-F2 C ) / ( DELY ( J ) +D ELY CJ ) ) 

F i A =V ( I , J + 1, K -1 ) / D EL / ( J ) • • 2 

p la=V( I , J , K-l ) • (1 . 0/ DFLY ( J ) A *2-1 .0/ DELY C-J )**2) 
c 1 C = V ( I , J + l, K-l ) / C E !_ Y ( J)' ; 2 

FI C=-Fi f 

f 1 = DE L Y ( J ) • Or L Y C J ) ( F 1 A +F 1 3-F1 C ) / ( DELY ( J ) +C EL Y ( J ) ) 

F= ( F2-f 1 ) / ( 2 . O l " 0 DLLZ) 

G=(W( I + 1 » J » K ) -W ( i — 1 1 J » K ) ) / ( 2 . 0 E0 5 i DELX) 
ri 1 = to( I , J+l ,K) /DE LY ( J) 

h2 = W( I, J,K)- ( l.OFO/QEL Y( J )+ 1. 06 O/DEL Y( J ) ) 

H2=W(I , J+l ,K)/CELY (J ) 

H - 2 . OF 0 • (H1-H2+H3) / ( D FL Y ( J ) +DE L Y ( J ) ) 

WTEE=-( A+B+r+O+RE I (=+=-G-H)) 
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to RITE ( 20, 904) « T E r 
C G.\T I NL - 

^LCJLATIC, jN THE TOP C'F THE CALCULATION CC GICN 
n- 2? I =1’ » L.f 1 

j~ 23 K = 2 , iv • 1 
J = £ 

J+l lEPMS HAVE BEEN ADJUSTED aITH UNIFORMITY ARGUMENTS 
r A lCULATION 0 F UT F R C V 3D MOMENTUM EQUATION 
A= <U( 1 + 1, J ,K ) -t -2-U { 1-1, J, K )•* ' 2) / ( 2. OEO* DELX) 

F 1 = 01 I , J,K ) V( I , J, KJ / (DELY ( J) ' 2) 

j 2 = U( I f Jf K ) V( : t j , K ) ( 1. OE O/DEL Y( J-l) = 2- 1 . OEO/ DE L Y ( J ) ~ 2 ) 

- ” = U l I t J-l , K ) ' V t I ,J-1,K)/“ELY(J — 1) ;: 2 
- l rL Y( J) * D E L Y ( J — 1 ) ■ ( 61+32-33) / < DELY ( J ) + CEL Y ( J-l ) ) 

C = (U( I* J, K + l) to ( I , J,K+1)-U( I, J , K- 1 ) W( I , J,K-1) ) /( 2. OEO? CELZ) 
D = < Phi < 1 + 1 , J,K) - PHI ( 1-1, J, K ) ) / < 2 .OEO* DELX ) 

E2A = V( I + 1, J,K) /DEL Y( J) • 2 

E2B = V ( I +1, J,K ) ( 1 .OEO/CELY( J- L)-'- 1 2-1.0E O/CEL Y( J) < *2) 

E 2C = V( I +1 , J-l , K) /DE LY{ J-l ) a " 2 

c 1- CEL Y(J) DELY(J-l) ( : 2A+E 2B-E2C ) /( DF LT ( J) +OELY ( J-l ) ) 

Ex^VI T-l , J, K J/GELY { J )■’ 2 

E IE -V( i-i , J,K) (1. JEC/DELY (J-l )•' 2-1 .0 EO / D ELY ( J-)- 2 1 
E 1 C = V I I - 1 * J- 1 1 K i / D F. L Y ( J - 1 ) •• • 2 

=1=DELY ( JJ 'DELY( J-l ) • ( eiA+=16-=lC)/ I DELY ( J ) +CELYIJ- 1) ) 

E = i E2- z 1) /( 2. 03 O c DELX) 

F=<n( I+lt J*K+l)-r.( I+i, J,K-1 J,K+1)+W(I-1, J, K-l) ) / 

IIh.OEO oel-x : DELZ) 

G 1= U( I , J , K ) /DEL Y( J ) 

G2=U< I , J, K )- ( 1 .3 EG /DELY ( J- 1 ) + 1 . 0 5 O/DEL Y < J ) ) 

G 3 - U ( I , J-l tK) /DE LY( J-l ) 

6 = 2. OEO ( G1-G^+G3)/(DFLY( J ) +CEL Y< J- 1 ) ) 

H = ( U ( I *J* K+l) -2. CEO U ( I» J,K )+U ( I, J,K-1 ) ) / DEL Z-" ■ • 2 
U”r E=- ( A + 6-+C + D + REI (E + F-G-H)) 

W- IT E C 19,963) UTEE 

F ALCL'LA TI ON r F WT =3^ 3D MOMENTUM EQUATION 

A = (U( I + lt J ,K ) ■ „ ( I + i, j, K )-U( 1-1, J ,K>’ 5 W( I - 1 , J,K) ) /( 2. OEO' DELX) 
ei=V(I , J,K) W( 1 , J ,K) /CELY ( J )r >2 

6 2= V ( I , J , K. } ■ to { I , J , K ) ( 1.0 C 0/DE LY( J-l) - 2- 1 . 0 K O/DE L Y ( J ) - 2 ) 

5^ = V< I, J-1*K)* W( I, J-1»K)/0=LY( J-l)**2 
a =: ELY( J-l )1 DE LY( Ji ( 61 + 62-63 )/ ( DELY ( J ) + CELY ( J-l ) ) 

C - ( to ( I , J , K I 1 ) •* 2-to ( I , J , K- I ) > ; 2) /( 2.0E0> DELZ) 

D = ( PH I ( I ,J,K+1)-?HI ( I , J , K- 1 ) ) / ( 2 .OEO' DELZ ) 
t = l U( I + 1 » J , K+ 1 ) - U( I +1 ♦ J » K — 1 ) -U ( I — 1 « J » K + l ) +U ( I— 1 » J , K-l ) ) / 

1 ( 4 . 0 c 0 DELX- DELZ) 

F 2 A = V ( 1 , J , K+ 1 ) / D E L Y ( J ) ’ 2 

F 26 = V ( ! , J ,K+1) ( 1. OE O/DEL Y( J-l ) ** 2-1 . OE 0/ DE LY ( J ) + - 2) 

F2C=V ( I, J-l, K+l J/CELY (j— 1) ' i 2 
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F?- Dr L i' (J ) *DELY { J — 1 ) ( F2A + r 2R- = 2C ) /( DE L Y ( J ) +DE L Y ( J- 1 ) ) 

p 1-=V( I , J » K-l ) /CELY< J ) • ; 2 

F IE = V< I , J t K-l) * ( 1. O/HE LY( J- 1) 2-1.0/DELY( J ) * " 2) 

F1C = V{ I, J-l, K-l ) / D t L Y ( J- 1 ) 

F 1 = OF L Y ( J ) 'DELY< J-l ) * ( F 1 A+ F 1 o-F 1 C ) / ( D^LY ( J) +OELY ( J -1 ) ) 

F= l F2-F 1) /( 2.0E0’ PBLZ ) 

G-(w(I+1,J,K)-a(I-1,J,K))/(2.J : 0 OFLX) 

H1 = W( I , J , K) /DE L Y ( J I 

r2=rt( I , J , K )~ ( 1.0=0/ CEL Y( J- 1) +1. OF 0/DFLY( J) ) 

H3=W( I , J-l ,K) /DELY ( J-l ) 

H= 2. OF 0 (H1-H2+H3) /( DEL Y( J) +QE LY( J-l ) ) 
w T E E = - ( A + B+C+D + e EI >( C +F-G-H) ) 

Rr ! TE (20,9041 WT EE 
CONTINUE 

CALCULATION ON THE CORNER 
CALCULATION OF UT FROM 3D MOMENTUM EQUATION 
CO 24 I =2 , LEI 
J-l 
K = 1 

J-l AND K-l T E - YS HAVE BEEN ADJUSTED 

A--( U( I + l» J ,K) • 2-U(I-l , J,Kh- 2) / (2.0E0- OELX ) 

p L=U( I ,J + 1,K )»V( I,J + 1,K)/(DELY( JK- - 2) 

Bz = U( I , J,K) + V( I , J,K) : ( 1 .OEO/UELY ( J|- *2-1 . 0 =0/ DE L Y ( J ) - 2 ) 

5 3- U( I , J+ 1 , K ) * V( I , J+l ,K) /DH L Y ( J ) : 2 

3=C?LY( J ) ~ OF L Y ( J ) •' ( 31 + B2-S3 ) / ( l, F LY ( J ) +OELY ( J ) ) 

C= ( U ( I , J , K + 1 ) w ( I » J , K + 1 ) - ; J ( I . J . K ) W ( I , J , K ) ) / ( 2 . 0 E 0 : D c L Z ) 
D=(PHI ( I + i , J,K J-PHI ( 1-1, J ,K ) )/ ( 2.0E0 - DEL X) 

E 2A=V( I + 1 , J+ 1 ,K) / D E L Y ( J ) * '• 2 

E2B = V'( I + 1 , J , K ) ( 1 .050/ DEL Y( J ) : ' 2- 1 . OE O/DE L Y( J ) i => 2 ) 
r 2 C = V ( I + 1 , J+ 1 , K ) / DEL Y ( J ) • 2 

1 2C=-r 2C 

E2 = DELY (J ) DEiY (J ) { E2A + E 2B-E 2C ) / ( DEL Y{ J ) +DELY{ J ) ) 

E1A = V(I-1 » J+ 1 » K ) / D E L Y ( J ) 4 *2 

E 13 = VI T-lt J,K ) • ( 1 . OEO/CEL Y( J) • 2- 1 . O c 0/ D E L Y ( J ) . •<- 2 ) 

E1C = V ( 1-1 ,J+1»K)/DELY(J)* ’• 2 
E 1C =- r 1C 

E1 = DELY ( J ) >DEL Y ( J ) ( E 1A + E 1B-F 1C ) / ( DFLY( J ) +OELY( J ) ) 

F = (52-:1)/(2.Q?0’ DELX ) 

F = ( W( !t 1 ,J , K+ 1 ) - J( I+1,J»K)-W(I-1 » J » K+l ) + W ( 1-1 ,J,K) )/ 
K4.0E0 DFLX'DSLZ) 

G 1 =U ( I , J+l ,K) /CELY( J) 

G2=U( I , J , K ) ( 1.0E0/DELY( J)+ 1.0E0/DELY( J) ) 

G3 = U( I , J + l ,K )/DELY ( J ) 

G = 2.0E O'- (G1-G2 + 33I /(DFLY( JI+DSLY ( J) ) 

H= ( U ( I,J,K+l)-2.0 r 0 U( I , J , K ) + U ( I ,J,K) ) / D E L Z 2 
UT = E=— (A+B+C+D+REI - (E+F— G— H)) 
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WRITE (19,903 ) 'J T 5 ■: 

r 

•: CALCULATION OF W T FROM 3D 30 T UM EQUATION 

A=(U(I+1»J,K)*W( I + 1 , J » K ) — U ( 1 — 1 , J . K ) ’ W ( 1-1, J ,N) )/ ( 2.0? U DEL X) 
= 1 = V( I , J+1,K) * W(I , J+ 1 , K) / DELY ( J ) : -t“2 

62-V( I » J, K ) w ( I, J,.<) ( l.OEO/DEL Y( J ) - 2- 1. 0E0/Dr.LY( J) 2) 

- 3 = V ( I , J + 1 » K ) M(I , J+l , Ki/ JELY ( J )•■ «'2 
s2=-H3 

b = DEL Y ( J ) D c L Y ( J ) ( 6 1+62-B 3 J / ( L t L Y ( J ) + 0 E L Y ( J ) ) 

C =( W< I , J ,K+1 ) «■' 2-rt (I t J , K) - 2) / (2 .OFC * OELZ ) 

D= ( PH ! ( I,J,K+1)-PHI(I,J,K)I / 1 2. OEO'DELZ) 

E = (U( I+i,J»K+l)-U( I+1»J,K)-U( I-1,J»K+1 )+U( 1-1 » J » K ) )/ 

1( 4.0 EO- DtLX DELZ) 

F2A=V(I,J + I,\ + l)/DtLY(J)»-i2 

r 2B = V ( I , J, K+l H ( 1 .OEO/DELY ( J) v , 2-1 .OEO/ DELY(J )- ^ 2) 
t ^C = V( I , J + 1, K, +1 ) /DE L Y( J) ' -2 
*-4iC = -F2 C 

P 2 = DEL Y( J ) ' OELY( J ) + ( F2 A+P2B-F2 C ) / ( OELY ( J ) + CELY ( J ) ) 

F 1A= V ( I ,J+1,K )/DELY( J) -2 

FI B=V( I , J ,K)~ ( 1 .0/ DELY ( J )1- - 2-1 .0/DELY( J I “ •* 2 I 
r 1 C = V ( I .J+1,K) /DD-LY( J) *•> 2 
F iC=-F 1C 

F1--D=LY( J) ' OELY ( J ) - ( FI £ + Fi 6-F1 C ) / ( DEL Y ( J J+PELYIJ ) ) 

F=( F2-F 1) /{ 2. 0: O' DELZ) 

G = (W( I + 1, J,K ) — to ( I - 1 , J , K ) ) / ( 2.0EQ : DELX) 

H 1 =W( I » J+ 1 ♦ K ) /OELY ( J ) 

H2=W(I,J»Ki‘( 1.0E0/DELY(J)+1.0E0/DtLY(J) ) 

H'. = W( I ,J+1 ,K)/DELY( J) 

H-2. OE 0 / ( H 1-H2 + H 3 ) /(DELY(J)+DELY(J) ) 

WTF c --( A+B+C+O+RE I s (E + F-^G-H) ) 

WRITE (20,904) wT EE 
2*» CONTINUE 
L 

f CALCULATION on a corner 

DC 2D I = 2, L c 1 
J^l 
K -NE 

C J-l AMD K+i TE=MS HAVE BEEN ADJUSTED 

C CALCULATION OF UT PS C-M 3D MOMENTUM EQUATION 

A = ( U( 1 + 1, J,K)**2-U(I-1, J,K) -• 2) / (2.0=0' DtLX) 

Bl-U( I »J + 1,K ) ‘ VU ,J + 1,K)/(0FLY( J )«■=•' 21 

b 2 = U ( I , J , K ) V( I , J, K) ( 1 .0 FO/ DEL Y ( J )**'*2- 1 .0 £0/ DEL Y ( J >•- 2) 

3 3 =U( I , J + 1 , K ) ' V ( I , J+l,^)/0cLY(J)-.'2 

o3*-B3 

rt =DE L Y ( J) s - DE L Y ( J ) ’ (B1+d2-B3 1/ ( DELY ( J ) + QELY ( J ) ) 

C = ( U ( I , J , K ) w ( I , J , K ) - U ( I , J , K- 1 ) ’ Mil . J.K-1) > /(2.0E0* DELZ) 
D=(PhI U+1,J,K )-PHI(!-l,J,K) )/( 2.0r0< , = LX) 
t 2A = V( 1 +1 , J+l , KJ /D C LY( J) * -2 
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E 2o = V ( I +1 t J i K ) *■ (1 .OEO/DELYU )" v 2 -1 .0 Eu/ DtLY(J )' >2) 

E 2 C = V ( I + 1,J + 1,K)/D5LY( JM* 2 
E2 C=-E2 C 

E2=0FLY< J) DELY( J) - ( r 2 «+ E 2S-E2 C ) /(DELY( J)+DELY(J ) 1 
E14=Vi I- I, J+1,K)/DELY< J M- 2 

Eifa = V( I -1 , J,K) M 1 .OEO/DELY ( J ) 2-1 .0 EO/ DELY ( J ) "2 ) 

r 1C=V( I-i, J+1,K) /DELY( J)-' 2 
~ 1 C=tE 1 c 

E1=DELY( J) ‘DELY( J) • ( El A + El S-El C )/ ( DELY ( J ) +0 ELY ( J ) ) 
c=( E2- = l)/( Z-OEQ-' CFLX) 

F=(w( 1+1, J,K)-W( 1+1, J,K-1)-W( I — 1 , J,K)+k( I-1,J,K-1))/ 

1( 4. OE 0; DELX' CELZ) 

G1=U ( I, J+1.K)/DELY( J) 

G2 = U( I ! J»K)’ (1 .0 EO/DELY (J )+l .OEO/DELY(J ) ) 

G 3 = U( I , J+l ,K) /CE LY( J) 

G= 2 .0 EO"' (G1-G2+G3)/(UELY( J )+DELY( J ) ) 

H = (U( I » J » K) -2 . OEO* U( I » J » K ) +U ( I » J » K-l ) ) / DELZ J - 2 
UTE3=-( A+B+C+D+REI- (E+F-G-H) ) 

FITE (19,903 ) UTEE 
C 

C CALCULATION OF WT F4QM 3D MOMENTUM EQUATION 

A = (U( 1+1 i JiK ) •+ W ( I + 1 » J * K ) — J ( I — 1 * J t K )•? W ( I - 1 » J » K ) )/( 2 • OEO * OF L X) 
3 1 = V( I ,J + 1»KM W(I ,J+1,K)/D=LY(JM*2 

22=V ( I ♦ J » K l v W ( I , J , K ) * ( l.OE O/DEL YUM '* 2- i . OE O/DE L Y ( J ) - ’ 2) 

£3 = V( I , J+l ,K)-W<;I , J+l » K )/ OELY ( J M*2 
B 3=-E 3 

B-DEl Y ( J ) ' DEL Y ( J ) : '(31 + B2-63)/(DELY( J)+DELY( J) ) 

C = ( W ( I i J»K) f • t 2 - W ( I , J . K-l )*••• 2 ) / (2.0 tO*Dc L Z ) 

D= (PHI ( I, J ,K )-3HI ( I , J,K-1) ) /( 2. OFO-DELZ) 

E=(U( 1+1, J,K)-U< I + 1,J,K— 1)— U(I — 1*J»K)+U( I-1,J,K-1) >/ 

1( 4. OE 0* L‘E LX DELZ) 

F2A=V( I , J + l, K J/DEL Y( J ) - > 2 

f 2 5 =V ( I » Jt KM < l.OEO/DFLYC J )* 2-1 . 0 tO/DEL Y ( J M > 2 ) 

F 2C = V ( ! ,J+1,K) /DELY(J) 2 
F2C=-F2C 

F ? = CELY < J) -DELY( J)> ( F2 A + F2 6-F2C )/ (DELY ( J ) + DELY ( J ) ) 

F 1A = V ( I , J + 1* K-l ) /PE L Y( J ) *»• 2 

C 1B=V( I , J,K-1 )■'•(! .O/DE I Y ( J ) •*' *2- 1 . 0/DEL Y ( J ) '- 2) 

F 1C « VC I » J+ 1 * K— 1 ) /Dt LY ( J) * 2 
F 1 C =— F 1C 

F1=DELY (J MDELY (J M { =1 A + F 1B-F1C ) / < DEL Y( J ) +DFL YU >) 

F - ( F 2— F 1) /( 2.0E0- ; DELZ) 

G=(W( 1 + 1, J,K)-k.( 1-1, J,K) )/( 2.0E0- ? DELX) 

HI = V« l I , J+l .K) /DELY ( J) 

H 2= W ( I t J * K ) -• ( 1.0EC/PtLY( J ) + 1. OE O/DE LY( J) ) 

H j =K ( I , J + 1,K )/DEL Y (J ) 

H =2. OE G J (HI-H2 + H3 ) / ( DFLY ( J ) +DELY ( J ) ) 
wTEE=-( A + B+C + D+REIME + F-G-H) ) 
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WR ITF ( 2C, 904) KTEE 
2 5 CONTINUE 
C 

C C :LCUL4'TIC"'I OF UT FROM 3D MOMENTUM EQUATION 

C CALCULATION ON A CORN p R 

C J+l AND K-I TERMS HAVE BEEN ADJUSTED 

D r ’ 26 1 = 2, LEI 
J = M E 
K = 1 

A = cue 1 + 1, J ,K )* 2-U ( I — 1 , J » K ) ; ; 2) / ( 2. OE O' C E LX ) 
81=U(I,J,K)-V(I,J,K)/( DELY ( J ) 2 J 

6 2=U( I , J ,K) V( I , J , K) - ( 1. OEO/DELY ( J-l ) - • 2~1 .UEO/CELY ( J) ' Z) 

B3 = U( I , J-1,K)< V ( I , J-1,K )/DELY( J - 1 ) “ '• 2 

b=DELY( J) v DELY( J-l )M B1+B2-B3 )/ ( CELY{ J )+DFLY(J-l ) ) 

C = { U ( I , J , K + 1 ) A W ( I , J,K + 1)-U( I to ( I ,J,K) ) / { 2 • 0 E 0 ' DELZ) 

D= (PHI ( 1+1, J,K )-PHI ( 1-1, J,K ) )/( 2.060*! DEL X) 

E2A = V( 1+1, J,K) /DELY( J)-^2 

E2d = V( 1+1, J,K)*( 1.0E0/DELY{ J- 1 ) ~ ' 2-1 . OE O/DE LY( J) "2) 

E 2 C = V ( 1+1 ,J-1,K)/DELY ( J-l ) ' •* 2 

E2=D5LY( J)-'DELY( J-l) '• { E 2A+ E 2B-E2 C ) / ( DEL Y ( J)+DELY { J-i ) ) 
E1A=V(I-1*J»K)/DELY(J )*-’ 2 

E16*V(I-1, J,K)M1.0E0/DELY(J-1 )*■ 2-1 .0 EO/ DE LY { J )<■> 2 ) 

E1C=V( 1-1, J-1,K) /DELY( J-1K*2 

El = utLY ( J ) ? DELY (J-l R ( E 1A+ E 1B-E 1C )/ ( DEL Y( J ) +DEL Y( J-l ) ) 
E=(E2- c l)/(2.0E0 i ‘DELX) 

R=(W(I+l,J,K+l)-to(I+l,J,K)-W(I-l,J,K+l)+w(I-l,J,K))/ 

1 K.OEO' DELX -'DELZ ) 

G 1 = U( I , J , K ) /DE L Y ( J ) 

GZ=U( I , J , K ). ( 1. OF O/DEL Y{ J-l ) + l. OEO/DELY( J) ) 

G3 = U( I , J-l ,K) /DELY (J-l ) 

G = 2.0E0> (G1-G2+G3) /(DELY{ J)+CLLY( J-l ) ) 

H= (U ( I , J , K + l ) -2 . 0E0 ; U ( I , J , K ) +U( I , J , K ) I /DE LZ* - 2 
UTE E=— ( A+ B + C+ D+R E I r - ( E + F-G-H ) ) 

WRITEt 19, 903) UTEE 

C 

C CALCULATION OF WT FROM 3D MOMENTUM EQUATION 

A=(LH 1 + 1, J,K )* to( 1 + 1, J,K)-U( 1-1 , J,K) W ( I - 1 , J,K) ) / (2 .OEO' DELX) 
6 1 = V ( I ,J»K)-to(I,J,K)/ DELY ( J 2 

B 2 = V( I , J , Kl ; W (I t- J , K) » ( 1 . 0 EO/DE L Y { J-l ) * 2 -1 . 0 EO/ D E LY f J ) - -Z ) 

B3=\A( I , J-1,K)*W ( I , J— 1 , K ) /DELY( J- l)** 2 

B=DELY( J-l )' DELY l J R ( B1+62-B3 ) / ( CELY( J )+DELY( J-l ) ) 

C = ( W ( I , J * K + 1 ) * v 2- to ( I , J , K) A 1 2 ) / ( 2 . OE 0* DE L Z ) 

D=(PHI( I,J,K + 1) — PHI(I,J-,K) )/( 2.0 EO* DELZ ) 

E = (U(I + 1, JiK + l)-U( 1+1 , J,K)-U(I-1,J,J< + 1)+U( I — 1 , J-, K ) ) / 

1 ( A . OE 0 ‘ DELX-' DELZ ) 

F2A=V( I , J ,K+1 )/DELY( J )**2 

FZE = V( I , J,K+1) * ( 1. OEO/DELY l J-l )*» 2-1 .OEO/DELY (J )+: 2 ) 

F 2 C = V ( I, J-l, K+l) /DEL Y( J-l) ^-2 
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F2= DEL Y ( J ) r E L Y ( J - 1 ) ( P-2A + F 25- F 20 ) / ( DE LV ( J) + Qt L Y ( J-l ) ) 

c 1 A =V ( I ,J,K.)/CeLY( J)"-‘2 

F 1 - = V ( I » J » K ) - ( 1 . U/DrLYU-l) 4 2-1 . C/DELY ( J )' •• 2 ) 

F1C = V( I ,J-1,K)/DELY( J-l)i' 2 

F1 = C5LY (J) * DE LY( J-l )*( FI A + FIE- F1C) / < DELY ( J }*DELY< J-l) ) 

F — ( F 2— F i) /( 2.0 : -0'DELZ) 

G=(w( 1 + 1, J,K )-,i( I-l,J,K))/(2.0E0f DELX) 

H1=W( I ,J,K)/DELY(J) 

H2=W( I , J, K ) ( 1. OF O/DEL Y< J- U+l. OEO/OELYt J) ) 

H3 = W( I , J-1,K)/DELY< J-l ) 

H = 2.0F0M H1-H2 + H3) /(DELY( J) +DELY ( J-l ) ) 

WT EE=-( A+5+C+D+R E I : - ( c + F-G-H ) ) 
wRITE (20,904) WTEE 
2c ^DNTINUE 
C 

C CALCULATION! _F UT FROM 3D MOMENTUM EQUATION 

C CALCULATION ON A CORNER 

C J+l A IN u K + l TERMS HAVE BEEN ADJUSTED 

DO 27 I =2, LEI 
J = ME 
K=NE 

A = ( U( I + 1, J ,K) r • 2-U( 1-1 , J, K) * - 2) /(2.0E0 DELX ) 

B1=D( I , J » K )'• V ( I rJ , K ) / ( DELY ( J-) - - 2) 

fc2-U(I , J , K ) - V(I , J,K)> ( 1 .OEO/ DELY (J-l ) • 2- 1 . CEC/ D EL Y ( J ) - 2 ) 

B 3=U( I » J- 1 »K ) *. VII , J- 1 , K ) /DE L Y( J- 1 ) 7 1 2 

9= DELY ( J ) ‘ DcLY ( J-l )■ ( B1+62-B3) / ( DEL Y( J) +DE L Y( J-l) J 

C = (U( I , J,K) ’ W( I , J ,K)-U( I, J,K-1 )N W (I , J,K-1 ) ) /( 2.0 EO- DEL Z ) 

D=(PHI ( I+1,J,K)-°HI< I-1,J,K) )/( 2. OEO 4 DELX) 

F2A=V( J+l»J»K)/DELY(J) a '2 

E2L = V( 1+1, J,K) ( 1. OEO/ LELY (J-l)* » 2-1.0 LO/C ELY (J) 1 * ^-2) 

E2C=V< 1+1. J-1,k)/DELY( J-l)**2 

E2 = DELY(J)i DELY (J-l )i (E2A+E2B-E2C)/( DEL Y(J )+OELY( J-l) ) 
tlA=V(I — 1,J,K)/DELY(J)**2 

E 1 b=V ( I - 1 , J , K ) ' ( l.OtO/OELYl J-1 )'m 2-1.0E0/DELY( J)v>-2) 

E 1C=V( I -1, J-l ,K) /DELY( J-l ) '2 

E1=DELY(J )"'DcLY( J-l) ■ ( E 1A + E 16-E 1C ) /( DE L Y ( J) +DELY ( J-l ) ) 
E=(E2- r l)/(2 ,OFD f DELX ) 

F = ( w( I + 1, J.K)— >,( I* 1 , J,K-1) - W( I -1 , J,K)+W ( 1-1 , J,K-1 ) )/ 

1(4. OEO DEl X ' DELZ ) 

G1=U( I , J, K )/ DE L Y ( J ) 

G 2 = U( I , J , K) ( 1. OEO/DF LY(J-1)+1.0E0/DELY(J) ) 

G3=U ( I , J-l,K)/r=LY(J-l) 

6=2.09 0*- (G1-G2 + G3 ) / ( DFLY(J ) +CELY (J-l ) ) 

H=(U( I , J ,K )-2. OEO* U( I , J,K) + U( I t J »K-1) ) /DELZ" *2 
U T E E = - ( A+B+C+D+RFI *• ( E+F-G-H ) ) 

WRITE (19,903) UTr c 
C 

C CALCULATION OF WT FROM 3D MOMENTUM EQUATION 
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c 

c 



c 



r 



A=(U( I + 1, J,K )*W( 1 + 1, J,K)-U( 1-1, J,K> *W( 1-1, J,K) )/(2.0E0t-DELX) 
31 = V( I , J,KH W( I , J, K) /DELY< J ) 

62= V( I , J ,K K w( I , J , K) ( 1. OEO/DE L Yf J- 1 ) * 2-1 . OEO/DE L Y ( J ) ♦ - 2 » 

63 = V( I, J-l ,K )"*W( I , J-l, K)/DELY( J-l)-='- 2 

B=DELY( J-l )> DELY( J) ♦ f B1+B2-B3) / ( DELY ( J ) + C ELY ( J-l ) ) 

C=(W( I, J,K 2- W( I , J,K-1) * : 2) /( 2. OEO*DFLZ) 

D = ( PHI ( ! , J ,K) -PHI ( I , J ,K-1 ) ) / ( 2 . 0 EO r ‘ DEL Z i 
E = ( U( I + 1,J,K)-U(I+1,J,K-1)-U(I-1 , JjKI+un-l t J»K-1) )/ 

1 (4.0E0- Dr LX DELZ ) 

F 2A = V( I , J, K) /DELY( J) >2 

F2B=V(I , J,K) M 1. OEO/DE LY( J- 1)“ • 2- 1. OEO/OE LY ( J )*■* 2) 

F2C=V ( i , J-l , K )/DELY (J-l )* u 2 

F2=OELY( J) >DE LY( J-l ) * ( F2 A> F2 B-F2 C ) / ( DEL Y ( J ) +DELY ( J-l ) ) 

F1A=V< I ,J,K-1)/DELY( J )*-> 2 

F iB=V ( I , J , K-l ) ’ ( 1.0/ DELY < J-l )* ' 2-1 .O/DEL Y ( J )v*2) 

F1C = V( 1 » J-1,K-1) /DELY( J-l>'‘-2 

F 1= DELY ( J ) E L Y ( J - 1 ) ’ ( F 1A + F IB- F 1C ) /{ DEL Y ( J ) +DE L Y ( J- 1 ) ) 

F = (F2-F1) / (2 . OEO* DELZ ) 

G= (W( 1 + 1, J ,K )-W( I- 1, J, K) ) / < 2. OEOt DELX) 

H1=W( I » J * K ) / DEL Y ( J ) 

H2=W(I , J,K)' (1. OEO/DE LY ( J-l )+1.0F0/ DELY (J) ) 

H3=W( I, J-l, K) /DEL Y( J-l) 

H=2 .0 EO* (H1-H2+H3 ) / ( DFLY ( J ) +DELY( J-l) ) 

WTf E=-( A + B+C + D+REI (E+F-G-H) ) 

W ^ IT E L20 , 904 ) WTEE 
2? CONTINUE 



CALCULATION AT ENDS OF THE REGION 
DP 28 J=1 , ME 
DO 28 K = 1 , NLE 

1 = 1 

UTEE=0. OEO 
WRITE( 19, 903) UTEE 
W T EE= 0 . C bO 
WR I TE ( 2 0 , 903 ) WTEE 
28 CONTINUE 



DO 29 J=1,ME 
DO 29 K=1,NE 
I =LE 

UTEE=0. OEO 
WRITE(19,903) UT EE 
WTFE =0. CEO 
WRITE ( 20, 903) WTEE 
29 CONTINUE 



IF(ITER.EO.O) DE LT = Df L T 2 
IF(ITEP.EQ.l) DELT = CELT 1 
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I F < IT EP .£0.0) Go TC 402 
I F(ICALC.EO.O) STD TG 402 
REKlNC 16 
REWIND 18 
RF 62(16) U 
R0 aD( 18 ) W 
4 0? CCNTINUF 
REWIND 19 
REWIND 20 



C -FAD VALUES GF UTEE AND WTR? AND COMPUTE NEW 

r 

DC 410 1=2, LEI 
CO 410 J = 2 » M E 1 
DC 410 K = 2 , N El 
R c AO ( 1 9 » 903 ) UTE? 

U( I,J,< )=U( It J ,K ) + UTEF* CELT 

R E AC(20 ,904 ) W7 E 5 

W( I , J ,K ) =W( I , J,K) + WTEF ; CELT 

410 COM INU E 

DC 411 1=2, LEI 
DC 411 J = 2 , M E 1 
K =1 

READ! 19 ,903) UTEE 
U(I,J,K. )=U(I,J,K ) +UT EF* DEL T 
E E AD ( 2 0 , 904 ) WTEE 
W ( 1 , J , K j = W ( I , J , K )+WTEE*DELT 

411 CCNTINUF 
r 

CC 412 1=2, LEI 
•DC 412 J = 2 , M E 1 
K=NE 

R FAD( 19, 903 ) U^E-E 
U<I,J,K)=U(I,J,iO+UTFE- DELT 
kEAD( 2 0 , 904) WTEE 
W( I,J,K) = W( I , J , K ) +WTf E'-‘ DEL T 
4 1 1. CCNTINUF 
C 

DC 413 1=2, LEI 
CC 413 K=2 , N2 1 
J = 1 

R EAD( 19 , 90 3 ) UTEE 
U( I , J, K) =U (I , J,K) + UT EE* DELT 
0 E AD( 20,904) WT^F 
/*( I,J,K)=W( I, J,K J+WTEE’’ CELT 
413 CCNTINUF 
C 



AND W 
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DO 414 1 = 2, LEI 
DC 414 K= 2 , N E 1 
J = ME 

PEAD(19,903) UTEE 

U ( I, J,K ) =U ( I , J , K ) +UT EF--DELT 

RFAD( 20,904.) WTE E 

W( I,J,K ) = W( I, J,K J+W'TEE- DELT 

414 CONTINUE 
C 

DC 415 1 = 2,1. El 
J = 1 
K = 1 

R E A D ( 19,903) UTEE 
U( I,J,K)=U( I, J,K)+UT EE- 4 CELT 
READ(20,904) WTEE 
W(I»J,K)=K(I»J,K) + WT EE-vDEL T 

415 CONTINUE 
C 

DC 416 1=2, LEI 
J = 1 
K = NE 

R E AD ( 1 9 » 903 ) UTEE 

U( I , J ,K) = U( I , J,K)+UTEE~* DELT 

READ(20,904) WTEE 

W ( I , J , K ) = W( I , J,K J+WTEE* DELT 

416 CONTINUE 
C 

GO 417 1=2, LEI 

J=F'B 

K= 1 

PEAD(19,903) UTEE 
U( I ,J,K)=U( I , J,K) + UTEE^ DELT 
P E AD ( 20 , 904 ) WTEE 
W(I,J,K) = W(I , J, K) +WT EE “DELT 

417 CONTINUE 
C 

DC 418 1=2, LEI 

J=ME 

K=NE 

Pt AD (19, 903) UTEE 
U( I,J,K) = U(I,J,K) + UTEE»DELT 
READ(20,904) WTEE 
W(I,J,K)=W(I , J,K)+WTFE S DELT 
A18 CONTINUE 
C 

DC 419 J=1,FE 
DO 419 K= 1 , N E 

1=1 
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419 



U(I, J,K ) = 1 .0 EO 
W ( I , J,K1 =0. OEO 

CONTINUE 






C 



c 

r 



C 

c 

c 



DC 420 J=1,MP 
30 420 K=1,NE 

U ( 1. E» J , < ) - 3 .0 E0 ’• U ( L r 1 , J , K ) — 3 . 0 E 0' U ( LE 2, J » K ) + J ( L5 3 » J , K ) 
w!(L = , J ,K) »3.0E 0+ W( LEI , J", K) -3.3S0 ' W<LE2» J, K ) +w (0*3, J, K ) 

420 CCi\T I N J fc 

DC 421 I = LP 1 ,L?2 
DC 421 K= NP 1, N° 2 

J = 1 

U( 1 , J , K ) =0. OE 0 
W( I,J,K 1=0. 0E0 

421 CCMINJE 
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CALCULATE V WITHIN REGION FROM 3D CONTINUITY 
C D 30 I =2 , LEI 
" 30 J = 1 , ME 1 
30 K = 2 , N E 1 

U C X2=( 0(1+1, J»K-1 ) +U( I+1,J + 1,K — 1 ) +U( I + 1 , J , K+ 1 )+U( I + 1,J + 1,K + 1) )/ 

14.0E0 

uEXl=(U(I— 1*J»K— 1 )+U( I— 1,J + 1,K-1)+U( I— 1,J»K+1)+U(I—1,J+1»K+1) )/ 



14. CEO 

0EX=(UEX2-UEXl)/( 2 . OEO-DE LX ) 



.;ZZ = (W( 1 + 1, J,K + l)+W(I+l,J + ltK+l)+W(I-l, J,K+1)+W( I- 1, J+l , K + l) I / 

14. OEO 

»Zl=(rt(I + l,J,K-l)+W<I + l,J+l,K-l)+W(I-l,J,K-lJ+rt(I-l,J+l,K-in/ 
14 .OEO 

WZE = < WZ2-WZ1) /( 2. OEO* CELZJ 
V ( 1 , J + 1,K )=V ( I, J , <<)-C&LY( J ) : (UE X + WZE ) 



CCr.TINUE 



CALCULATE V ON A SIDE FROM 3-C CONTINUITY EQUATION 
K — 1 TERMS HAVE 5cEN AC JUST ED 
DC 31 I =2, LEI 
CO 31 J = 1 , ME1 
K = 1 

0 E X 2= ( U ( I + 1,J,K) + U( I + ltJ+ltK) +U(I + 1, J,K+1)+U(I + 1, J + l, K + l) )/ 

14. OEO 

UEX1 = ( U( I“1,J,K)+U(I~1,J+1,K)+U ( 1-1 , J,K + 1 )+U( 1-1 , J + l, K + l ) )/ 

14. OEO 

UEX=(UEX2-UEX1)/(2.0E0 DELX ) 

Wi.? = ( n( I + 1 , J , K+l) +WU + 1 , J+l ,K+1 )+M 1-1 , J , K+1)+W(I-1 , J+l, K + l ) )/ 
14. OEO 

WZ1 = ( w( 1 + 1 , J ,K) + W( 1 + 1 , J+l , lO+rt ( I -1 , J,K) +W ( 1-1, J + l, K ) )/ 

1 4 . C E 0 
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non o o oo oo on o on 



WZE = < WZ2-WZ1) /( 2. 0E0' DELZ) 

V{ I , J + lt K)=V( It J,K )-DELY( J ) MUEX+WZE) 

31 CONTINUE 

CALCULATE V GN A SIDE FROM 3-D CONTINUITY EQUATION 
K+l TERMS AH Vt BEEN ADJUSTED 
DO 32 1=2, LEI 
CC 32 J = 1 , MF 1 
K = NE 

UEX2=( U( I + i, J,K- 1 ) +U( 1 + 1, J+ 1,K- 1) +U( 1 + 1 1 J« K)+U(I +1 t J+l f K)) / 
14. CEO 

UF Xl = ( U( 1-1, J, K-l) +U( 1-1 , J+l ,R-1 )+U( 1-1 , J , K)+U( 1-1 , J+l , K) )/ 
14.0E0 

UEX = (UEX2-UEXl)/(2 .0 EO 1 ' DELX ) 

W Z 2 = ( X ( 1+1, J,K)+W( 1+1 , J + l ,K)+W( 1-1, J,K) +W(I-1, J+l ,K) )/ 

1 4 . 0 EG 

WZ1=( W( 1 + 1 , J, K-l )+W(I+l , J+l ,<-l )+W(I-l , J, K-l )+w ( 1-1, J + l, K-l ) )/ 
14.0E0 

WZE=(WZ2-WZ1 ) / ( 2 . 0 EO ' DEL Z ) 

V( I , J+ 1 , K ) =V( I , J,K)-DELY( J) * (UEX+WZE) 

32 CONTINUE 

SET V EQUAL TO ZERO IN THE PLANE OF THE PLATE 
DC 34 1=1, LE 
DC 34 K = 1 , N E 
V( I,1,K) = 0.0E0 

34 CONTINUE 

SET UPSTREAM V EQUAL TO ZERO 
DC 33 J =1 , MF 
DO 33 K =1 , NE 
V< 1, J,K ) = 0 . OE 0 

35 CONTINUE 

EXTRAPOLATE DOWNSTREAM V FIELD 
DC 36 J=I ,ME 
DC 36 K = 1 , NE 

V(LE,J,K) = 3.0E0*V(LE1, J,K)-3.0E0* V(LE2, J,K)+V(LE3, J,K) 

36 CONTINUE 

WRITE (6, 6 12) 

F 12 <=CFMAT( 1H0,20X, • INTO CALCULATION POISSON') 

200 CONTINUE 

CALCULATE PRESSURE FIELD FROM POISSONS EQUATION 
ERR = 0. OE 0 

CALCULATION OF PRESSURF FROM 3D POISSON'S EQUATION 
CALCULATION WITHIN REGION 
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Of: 37 I =2, LEI 
D r 37 J = 2, ME 1 
CC 37 K=2 » NE1 

A 1= ( P HI (I+1,J,K)+PHI(I— 1,J,K)) /DbLX*- -2 
A2=(PHI (I, J,K+1)+PHI (I, J,K-1) )/DELZ7>2 

A 3 =2 . 0 E 0" (PHI (I , J+l,K)/DELY(J)+PHItI,J-l,K)/DELY( J-l)}/ 

1( L'ElY( J )+DELY( J-l) ) 

A = A 1+A2+A3 

3 = (U(I + l,J,K)**2-2.0E0 4 U(I,J,K)^2-rU( I- 1 , J , K )>- 2 ) / DEL X - <-2 
C1=V( I t J + 1»K 2/DELY( J) 

C2 = (V( I ,J,K)* :4 2 )‘( 1.0 EO/ DELY (J- l)+i. O/DEL V(J) ) 

C 3 = V( I »J-1}K)'"2/DELY(J-1< 

C = 2.0E0-MCi-C2+C3) / (DELY( J )+DEL Y( J-l) ) 

D = ( W(I , J,K*1) *'2-2 .OEO «W( I , J ,K )* - 2+W( I, J , K- 1 ) a -- 2 ) / DELZ :j - s 2 
E 2A = U( I + lfJ+l»K)' V(I + lrJ+i»K)/DELY(J) < 2- 

F 2 E=U ( I + 1 » J, K ) • V ( I + 1 » J , K ) M 1 . 0 E 0 / D E L Y ( J - 1 ) 5 *’ 2- 1 . OF O/OE L Y ( J ) -- - 2 ) 
E 2C=U ( I + 1 » J-l » K) • V ( I +1 ,J-1 , K-)/DELY( J-l) ^^2 
E 2= DEL Y ( J ) * DEL Y ( J- 1 ) ( E 2A + E 2B-E 2 C ) / ( DE L Y ( J) +DELY ( J-l) ) 

51A=J( 1-1, J+1,K)«V(I-1,J+1,K)/DELY( J)**2 

E1B=U( 1-1 , J,K)< V( 1-1 , J,K) M 1.0E0/CELY( J-l )**2-l . 0 EO / DEL Y ( J ) 2 ) 

F 1C=U ( 1-1, J- 1,K ) +V( I- l t J-l, K) /DELY( J-l) 2 

E1 = DELY ( J )*DELY (J-l) - ( E1A+E 1B-E 1C ) / ( DEL Y ( J ) +0 EL Y( J- 1 ) ) 

B = ( E 2- E 1 ) / ( 2. OE 0* DE LX ) 

F2A=V ( I»J + 1,K + 1)' T W( I,J + 1»K+1) /DELY(J) £ a 2 

P2B = V( I , J, K+l )*W ( I , J,K + 1 ) ( 1 .0 50/ DELY ( J - 1 )*!••' 2-1.0£0/DELY(J ) •' • 2) 

p 2C=V ( I , J-1»K+1)¥W( I » J-1,K+1) /DELY< J-l)**2 

F2 = DELY ( J ) :) DELY ( J - 1 ) + ( F2A + F 25- F2C)/(DELY(J)+DELY(J-1) ) 

F1A = V( I , J+1,K-1 ) v W( I , J+l , K-l )/DELY( J)** 2 

F1B=V( I ,J,K-1)‘ U( I , J ,K-1) M 1.0E0/DELY( J-l)*-* 2-1 . OFO/ DELY ( J)- ■: 2 ) 
F1C = V( I , J-I , K-l ) ■*■ W ( I, J-l, K-l)/ DELY (J-l)** 2 
F i = Dt L Y( J ) ! 0 E L Y ( J-l)- (F1A+F13-F1C)/(DELY(J) +DELY (J-l ) ) 

F= ( F2-F 1 )/ ( 2.0F0-DFLZ ) 

G2 = U( I + i, J,K+1 ) * W( 1+1 , J,K+1 i+U< 1-1, J,K-1 )>W ( 1-1, J, K-l) 

G1=U( I-1,J,K+1)* W( I-i, J,K+1)+U( I+1,J,K-H*,W(I«-1, J,K— 1) 
G=(G2-G1)/ ( 4 .OEO ~ DELX*DEL Z ) 

H = (2. 0E0/DELX*>-2+2 .OEO/DELX’ 4 -2 + 2 .OEO/ ( DELY ( J J + DEL Y ( J-l ) ) ) 

A p HI= ( l.OEO-GMEGA J---PHI ( I ,J,K)+CMEGAr ( A+ B + C + C+2 . 0 EO- 1 (E+F + G) )/H 
RcSIC-ABS ( APHI-PHI ( I, J , K) ) 

FCC=m£X1(E d R,RESIC) 

PHKI, J,K)=APHI 
3 7 CONTINUE 
C 

C CALCULATION OF PRESSURE F»3M 3D POISSON • S EQUATION 

C CALCULATION ON A SIDE 

C K-l TERMS HAVE BEEN ADJUSTED 

DO 3# 1=2, LEI 
DC 38 J =2 , ME1 
K = 1 
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c 

c 

c 

c 



A1=(PHI ( 1+1, J,K)+PHI( I~l» J,K))/DELX**2 
A 2= (PHI ( I « J f K + 1) + PH I (I fJ«K) ) /DELZ**2 

43=2. OEO* '(PHI ( I, J+ltK) /DELY( J)+PHI (I , J-l ,K) /DELY( J-l) )/ 

1(CELY(J ) + DEL Y( J-l) ) 

A =A 1+ A2 + A3 

b = (U( 1 + 1, J,K)* t 2-2.0E0* U( I , J,K)**2+U(I-1 , J , K) ** 2 ) / DELXM 2 
C1=V( I , J+l ,Km 2/DELY ( J ) 

C 2 = ( V ( I , J, K) **2) *(1.0E0/DELY( J-l 1+1.0/ DELY (J) ) 

C3=V( I , J-1,K)**'2/DELY( J-l) 

C=2.0E0v(Cl-C2+C3)/ (DELY( J )+DELY( J-l) ) 

D=( W( I , J,K+1) **2-2. OEO*W(I , J,K)**2+W(It J , K )**2 ) / DELZ**2 
E2A= J ( I+ltJ+ltKJ’vV (I+1»J+1»K) /DELY(J)£‘2 

62S=U(I+1»J»K)*-V(I+1»J»K) I M1 .OEO/DELY( J-l )**2-l. OEO/DELY( J )**2) 

c 2C=U( I+1,J-1,K)*V(I+1,J-1,K)/DELY(J-1)**2 

E2 = DELY (J ) f "DELY (J-l )* ( E2A+E 2B-E 2C ) / ( DEL Y ( J ) +DE LY( J - 1 ) ) 

E 1A=U( 1-1 t J+l ,K) *V( I— 1 , J+l , K)/DELY ( J)**2 

£ 1B=U ( I-1,»J,K)*V( 1-1, J,K)*( 1.0E0/DELY( J- 1 ) * * 2-1 . 0 EO/ DELY ( J )** 2 ) 
El C=U ( 1-1 , J-l ,K)*V( 1-1 , J-l t K)/DELY(J-l)** 2 
E 1=DE L Y ( J) OE LY( J-l) * ( E 1 A+ E 16-E 1 C ) / ( DE L Y ( J ) +DELY ( J -1 ) ) 
E=(E2-E1)/(2.0E0*DELX) 

F2A=V ( I ,J+.1,K+1) *W(I t J+1,K+1 )/DELY(J )**'2 

F 2B = V( I ,J,K+1)*W(I ,J,K+1)M1.0EO/CELY( J-l )^2-l . 0 EO/ DELY ( J )** 2 ) 

F2C = V ( I , J-1,K+1)*W( I, J-t,K+l) /DELY( J-l)**2 

F2=DELY( J)*DELY( J-l )* (F2A+F2B-F2C)/ (DELY (J )+DELY (J-l ) ) 

F 1A=V ( IfJ + l»K) + W( I »J+1»K)/DELY( J ) **2 

F1B=V( I , J,K)*W(I t J,K)5 ( 1. 0 EO/ DELY ( J- 1 )** 2- 1.0E O/DEL Y( J )**2) 

F1C = V( I ,J-1,K)*W(I t J-l »K)/DELY( J-l)* *2 

F1=DEL Y (J)*DELY( J-DMF1A + F 1B-F 1C) / ( DEL Y ( J ) +DE LY ( J-l ) ) 

F = (F2-F1)/ (2.0E0J-DELZ) 

G2=U( I+l,J,K+l)lW(I+l,J,K+l)+U(I-l,J,K)*W(I-l,J,K) 

G1=U( 1-1, J t K+l)tW( 1-1, J,K+1)+U( I+1,J,K)*W(I+1, J,K) 

G = (G2-G1) /(4.0E0*0ELX*DELZ) 

H= (2.0 EO/D EL X*» 2 + 2. 0E0/DELX* 4 2+2. 0E0/( DELY (J)* DELY ( J-l ) ) ) 

APHI= ( 1 .OEG-OMEGA)*PHI ( I , J , K ) +0 MEGA* ( A+B+C+D+2. OE 0+(E+F+G) ) /H 
RESID=ABS( APHI-PHI ( I * J t~K) ) 

EPR=MAX1( ERR, RES ID) 

PHI M » J » K ) =APHi 
30 CCNTINUE 



CALCULATION OF PRESSURE FROM 3D POISSON'S EQUATION 

CALCULATION ON A SIDE 

K+l TERMS HAVE BEEN ADJUSTED 



DO 39 I =2, LEI 
DO 39 J=2,ME1 
K = NE 

A1=(PHI (I + ltJfK) + PHI ( I-l»J*K))/DELX+*'2 
A2= (PHI ( I,.J,K)+PHI( I, J,K-1) ) /DELZ^^2 

A3=2. OEO* ( PHI (I» J+1,K) /DELY ( J) +PHI ( I , J-1,K )/DELY ( J-l) )/ 
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1 ( LELY ( J )+CELY ( J-l ) ) 

A =A 1 + A 2+A3 

B= (LH I + 1, J,K )< «■ 2-2.0t0*U( I , J4 ,<) * ^2+U( I - 1 , J ,K)*+2) /DELX* *2 
C1=V( I , J+l ,K)M 2/DELY ( J ) 

C2 = (V< I ,J,K) **£)"*( 1.0E0/DELY(J-1)+1.0/DELY(J) ) 

C3 = V ( I , J-l, K)*« 2/DELY (J-l ) 

C =2. OE 0* ( C1-C2+C3 ) / ( CELY ( J ) + DELY ( J-l ) ) 

D= ( U ( I , J,K) **2-2.UE0‘ W( I , J ,K)$ < 2+WU , J , K- 1 ) ** 2) / DELZ*/) 2 
E2A=J ( 1+1 , J+l t K) YV( 1 + 1 , J+l, K)/DELY( J )*>2 

E 2B-U( 1+1, J,K)1V( 1+1 , J,K) U 1.0 EG/ CELY (J-l K*2-l .OEO/DELY (J )*-*2) 

E2C = U( 1 + 1, J-l,K)rV( I + 1 , J- 1, K) /DE L Y( J-l ) M 2 

E2 = DEL Y ( J )" : DE LY ( J-l ) * ( E2A+ E 2B- E 2C ) / ( DEL Y ( J ) +DEL Y( J-l)) 

E 1A=U( 1-1, J+1,K) *V(I-1,J+1 ,K)/DELY( J)v*2 

E1B=J ( 1-1, J,K)*\T( 1-1, J,K)*( 1 . OEO/DELY( J- 1)**2-1. OE O/OE L Y ( J ) *■* 2 ) 
E 1C=U( I -1 , J-l ,K) *V ( 1-1, J-l , K )/DELY( J-l )*-2 
E1 =DElY(J) ! * 0ELY( J-l)t (E 1A+E 1S-E 1C ) / ( DE L Y ( J ) +DE LY ( J-l) ) 

E= ( E2-E1 )/ (2 .OEO^GELX ) 

F2A=V(I , J+1,K)«W( I ,J+1 ,K)/DELY( J)**2 

F2b=V( I ,J,KH Wl I , J,lO*( 1. OEO/DELY ( J-1)M 2-1 . OEO/ DE LY ( J ) **2 ) 

F2 C = V ( I , J-l ,K ) ( I , J-l ,K)/QELY ( J-l)**2 

F2=DELY( J)“ 3ELY( J-i) b ( F2 A+ F 26-F2 C ) / ( DEL Y ( J ) +DELY ( J -L ) ) 

FiA=V( I ,J + 1,K-1)*W( I, J + 1,K-1)/DELY( J)**2 

Fie=V( I , J» K-l ) * w ( I , J,K-i )* ( 1 .OEO/DELY (J-l )*4 2-1 .OEO/DEL Y( J ) * ^ 2 ) 
F 1C =V( I * J- i» K- 1) F W ( I , J-1,K-1)/DELY( J-l)**2 
Fl = DcLY (J)*DELY(J-1 )< ( F 1A + F 1B-F 1C ) / ( DEL Y ( J ) +DE LY ( J-l ) ) 
F=(F2-F1)/(2.02C*DELZ) 

G2=U( 1 + 1, J ,K)« W( 1 + 1, J ,K) + U( 1-1 , J,K-1) *K( 1-1 , J,K-1 ) 

G1=U( I — 1 , J,K)*M 1-1, J,K)+U( 1 + 1, J,K-1)*W( I+1,J,K-1) 

G = (G2-G1) / (4.0E0i-DELX*'CELZ) 

H= (2.0EO/CELX$-A2 + 2.0£0/DELX+*2+2.0EO/( DELY( J)*-OELY( J-l) ) ) 
APHI=(1.0 EO— GMBG A ) + PH I ( I , J , K ) +OMEGA* ( A +B +C+D+2 .0 EO* ( E+F+G) )/H 
RESID=ABS( A PH I -PH I (I , J ,K) ) 

ERR=MAX KERR, RES I D) 

PHI (I , J ,K) =APHI 
39 CONTINUE 
C 

C C/-LCULATI CN OF PRESSURE FROM 3D POISSON'S EQUATION 

C CALCULATION ON TOP 

C J+l TERMS HAVE 3EEN ADJUSTED 

DO 40 I =2, LEI 
DC 40 K = 2» NE 1 
J =ME 

A1 = (PHI ( I + 1,J,K) + PHI ( 1-1, J,K) ) /DELX**2 
A2= (P H I ( I,J,K+1)+PHI(I,J,K-1))/DELZ**2 

A3=2.0EG> ( PHI ( I , J . KJ/DELY ( J )+PH I ( I , J-l ,K )/DELY( J- 1) )/ 

1( 0ELY( J)+DELY( J-l) ) 

A = Al+A 2+A3 

B-( U( 1 + 1 , J ,K) t -* 2-2 .GEO* U( I , J,K)$*2+U ( 1-1 , J , K)**2)/ DELX**2 
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c 

c 

c 

c 



C 1 = V( I , J,K)**2/DELY( J ) 

C2=(V( I , J,K) A *-2)*=( 1.0E0/DELY(J-1)+1.0/DELY(J) ) 

C3=V ( I» J-1,K )*%-2/DELY ( J — 1 ) 

C=2.0EOi (C1-C2+C3 ) / ( DELY( J ) +DELY ( J-l ) ) 

D=( W( I » J,K + l)*-4 2-2.0E0* W( I , J,K)**2+W( I , J , K-l)**2) /DELZ**2 
E2 A=U ( 1+1, J,K)*V< 1+1, J,K)/DELY( J )**2 

E 2B-U( I+1,J,K)*V(I+1,J,K)* ( 1.0E0/DELY( J-l )**2-l .0 E0/ DELY ( J )**2) 

E2C=U ( I + 1,.J-1,K)*V(I + 1,J-1,K)/DELY(J-1)**2 

E2=DELY ( J)*DELY( J-l )*■ ( E2 A + E2B-E2 C )/ ( DEL Y ( J )+DELY( J-l)) 

E1A = U( 1-1, J,K)» V( 1-1, J ,K) /DELY( J)**2 

E1B=U ( I — 1 , J,K)*V ( 1-1, J,K)*( 1.0E0/DELY( J-l) **2-1. OE O/DEL V(J).)*2) 

E 1C=U( 1-1 , J-l ,K)*V( 1-1 , J-l ,K)/DELY ( J-l )**2 

E1=DELY( J )*DELY( J- 1)* ( E1A + E 1B-E 1C ) / ( DE LY ( J ) +DELY ( J-l ) ) 

F=(E2-E1)/(2.0E0*DELX) 

F2A = V( I » J ,K>1 ) ♦ W( I , J,K+1) /DELY { J)**2 

F2B=V( I ,J,K+1 ) * a) ( I» J,K+1)*( 1. OE O/DEL Y( J- 1 ) ** 2-1 . OE O/DE LY ( J ) ** 2 ) 

F2C=V ( I , J-l ,K+1) *W( I , J-1,K+1 )/ DELY (J-l )**2 

F2=DELY( J)*DELY( J-l)* (F2A+F2B-F2C) /( DELY ( J J+DELY (J-l) ) 

F1A=V( I , J,K-1)*W( I, J,K-1)/DELY( J )**2 

F 1 B = V ( I ,J,K-1)N(I ,J,K-1)» (1 .OEO/ DELY (J-l )$*2-1.0E0/DEL Y( J )**2) 

F 1C=V ( I ,J-1,K-1)*W( I , J-1,K-1)/DELY( J-l)** 2 

F1 = DELY (J )'*DELY (J-l )* ( F1A + F1B-F1C)/(DELY(J )+DELY( J- 1) ) 

F=(F2-F1) /(2.0E0*DELZ) 

G2=U( 1+1, J,K+1) *W( 1+1, J,K+1)+U( 1-1, J,K-1)*W(I-1, J , K-l ) 

G1=U( 1-1, J,K+1) *W( 1-1, J,K+1)+U( I+1,J,K-1 )*W(I+1, J,K-1) 

G=( G2-G 1) /( 4. OE 0* DELX* DELZ) 

H=(2.0EO/DELX**2+2.0EO/DELX**2+2.0EO/(QELY( J)*DELY( J^l) ) ) 

APHI = ( l.OEO-OMEGA)* PHI ( I , J » K ) +OMEGA* ( A+B+C+D+2. OEO* ( E+F + G) )/H 
RESID=A8S( APHI-PHI ( I ,~J ,K) ) 

ERR=MAX1 ( ERR, RES ID) 

PHI { I , J ,K ) =A PH I 
40 CONTINUE 



CALCULATION OF PRESSURE FROM 3D POISSON'S EQUATION 

CALCULATION ON THE BOTTOM 

J-l TERMS AHVE BEEN ADJUSTED 

DC 41 1=2, LEI 

DO 41 K = 2 » NE 1 

J = 1 

A1=(PHI(I+1,J,K)+PHI(I-1,J,K)) /DELX**2 
A2=(PHI ( I, J,K+1 )+PHI ( I, J,K-1))/DELZ**2 ,,, 

A 3=2. OE Ot( PHI ( I , J+l ,K) /DELY ( J) + PHI H * J+l , K ) /DELY ( J ) \J 
1 ( DELY ( J ) + DElY( J ) ) 



A=A1+A2+A3 

B = (U( I + l« J,K)**2.-2. OEO* U( I , J,K)*+2+U( I~1 , J,K)**2 )/DE1_X*»2 
C1=V( I,J+l,K)*+2/ DELY ( J ) 

C2=( V( I , J, K) 4*2) * ( 1 . 0 EO-/DELY ( J ) +1 .0/ DEL Y ( J ) ) 

C3=V( I»J+1»K)*^2/DELY(J) 
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C a =— C 3 

C=2.0E0> (C1-C2+-C3 )/(DELY( J)+DELY( JH 

D = ( lnl( I , J,K+1 J * '2-2 . GEO* W ( I , J , K ) **2+W ( I , J , K- 1 )**2 ) / DEL Z*-* 2 
E 2A=U ( I+1,J+1,K) S>V( I + 1,J+1,K)/DELY( J)**2 

E2 E=U< I+l,J t K)-:V( 1 + 1, J,K)*( l.OEO/DELYU )** 2-1 . DE O/OEL Y( J)**2) 
E2C = U( I+1,J+1,K)*V(I+1,J+1,K)/DELY(J)**2 

E2C=-E2C 

E2 = DE LY ( J) *-DELY(J)ME2A + E2B-E2C)/(DELY< J )+DELY(J) ) 

E1A = U( I - 1 » J+ 1 , K ) <t V( I-l»J+l»K)/DELY(J)i <k 2 

tlB=U( 1-1, J,K)> V( I -It J, K) fc ( 1.0t0/DtLY( J ) ** 2-1 . OE O/DEL Y ( J)f^2) 
E1C=U< 1-1 , J+l ,K)« V( 1-1 ,J+1 ,K)/DELY(J)**-2 
E 1C=-E 1C 

El = DELY (J )*DELY ( J )4 ( El A+E 1B-E 1C ) / ( DELY( J J+DEL Y( J) ) 

E = ( E2- E 1) / ( 2 . OEOA DE LX ) 

F2A=V( I f J + l t K+l)tW( It J + ltK+l) /DELY(J)«*2 

F2B=V( I , J,K+1)*W( I , J,K+1) v ( 1.0E0/DELY( J )**2-l . OEO/ DELY( J )$*2) 
F2C=V( I , J+1,K+1) *W(I t J + ltK+1) /DELY( J)**2 

F2C=-F2C 

F2=DELY(J)*DELY( JIM F2A+F2B-F2C) / (DELY (J)+CELY(J ) ) 

F1A=V( I ,J*1,K-1)*W(I » J + 1,K-1J /DELY( J)*‘ 2 

FI 6= V ( I , J , K-l )*WII,J,K-l)t( l.OEO/DELYC J )*» 2-1 . OE O/DEL Y( J)**2) 
F1C = V( I ,J+1,K-1)*W(I » J+l, K-l) /DELY (J)**2 
F1C=-F 1C 

F1=DELY ( J) *■ DELY ( J)* ( F 1 A + Fl B-Fl C ) / ( DEL Y ( J ) +DELY ( J ) ) 

F = (F2-F1)/(2.0E0*DELZ) 

L>2=U ( I + lt JtK + l)«W( I + lt J,K*1)+U< I-lt JtK-1)* fc(I-l, Jt K-l) 

G1 = U( I -1, J,K+1) ”W( 1-1 , JtK + 1 ) +U{ 1+1, J,K-1 )*W (1 + 1, Jt K-l) 

G= ( G2— G 1 ) / ( 4. OE 0**D EL X*'DE LZ ) 

H= (2 .0E0/CELX**2+2 .0 EO / DELX* v 2+ 2 . OEO/ ( DEL Y( J ) ? DEL Y ( J ) ) ) 

APHI = ( 1 . OE G-GME3A ) ■*' PHI (I ,J,KJ+ ONEGA* (A+B + C + C+2 .0 EOM E+F+G) 1/ H 
RES ID-AbS ( APHI-PHK It J t K) ) 

ERR=MAX1(ERR, RES ID) 

PH I ( I , J ,K) =A PH I 
4L COKT I NU E 
C 

C CALCULATION OF PRESSURE FRGM 3D POISSON'S EQUATION 

C CALCULATION ON A CORNER. 

C J-l AND K-l TERNS HAVE BEEN ADJUSTED 

DO 42 1=2, LEI 
J = 1 
K = 1 

Al = (PHI ( 1 + 1, J,K)+PHI( 1-1, J,K)) /DELX**2 
A2 = (PHI (I , J, K+l ) + PH I ( I , J , K ) T/ DEL Z* s 2 

A 3=2. OE G* C PHI ( I , J+ 1,K) /DELY( J) + PHI ( I , J + 1 , K ) /DE LY ( J ) ) / 

1 ( D ELY ( J ) +DELY ( J ) ) 

A = A 1+ A2 + A3 

B = (U( I + 1, J , K ) 2-2.0E0 t U( I , J , K. ) *■*- 2+U( I - 1 , J , K) *+ 2 ) / DELX+ »2 

Cl = V ( I , J + l , K )*• f-2 / DELY ( J ) 
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c 

c 

c 

c 



C2= ( V ( I t JtK)**2 )*( l.OEO/D EL Y( J ) +1. O/DEL Y( J ) ) 

C3=V(I , J + l r K)#*2/DELY< J) 

C 2=— C3 

C =2 .OEO* (C1-C2+C3 )/ (DELY(J)+DELY(J ) ) 

D = ( W( I , J,K+1)**2-2.0E0* W(I , J ,K) **2+W < I , J , K )**2 )/DELZ*+2 
E2A=U ( I + 1,.J+1,K)*V< I + l» J+1,K) /DELYC jm 2 

£2B=U( I +1 ».J, K) ?V(I+i,J,KI>ll .0 EO/ DEL Y( J )** 2-1 .OEO/ DELY( J )^^2 ) 
E2C = U( I + 1,.J+1,K)*V(J + 1,J+1 ,K)/DELY( J)**2 

E2C=- E2C 

E2=DELY(J)’ t, DELY(J)»’(E2A + E2B-E2C)/(DELY(J)+DELY(J) ) 

E iA=U ( 1-1, J+1,K)£V( I-1,J+1,K) /DEIY( J) M2 

E16=U( I-l,.J f K) + V(I-1, J f K)*< 1.0EC/DELY(J )**2-1.0E0/DELY( J)**2) 

E 1C = U( 1-1 , J+l t K) £ V( 1-1 , J+l , K) / DELY ( J)**2 

ElO-EIC 

E1 = DELY (J )*DELY ( J )* ( E1A + E1B-E1C )/( DELY( J ) +DEL Y( J ) ) 

E=(E2-E1) /( 2.0EO+DELX) 

F2A=V( I, J + 1,K+1HW( I» J + l t K+l)/DtLY( J)4<*2 

F2 B-V ( I ,J,K+1)*W(I,J,K+1)*(1.GE0/CELY(J )**2-1.0E0/DELY(J )+*2) 

F 2C— V ( I t J+l»K+l)*W(I ,J+1,K+1)/DELY( J)**2 

F2C=-F2C 

F2=DELY( J)*DELY( JM- (F2A+F2B-F2C)/ ( DELY ( J ) +DELY ( J ) ) 

F 1A-V ( I t J.+ lf K)»W( I « J + 1,K) /DELY( J)**2 

F1B=V( I ,J,K)*W( I, J,K)*(1.0E0/DELY(J)**2-1.0E0/DELY( J )**2) 

F 1C = V( I » J+l »K) i W( I ,J+1 ,K)/DELY( Jl**2 
F 1C=-F 1C 

Pl = OELY (J )*DELY ( J )*< F1A + F1B-F1C )/ (DELY( J )+DELY(U ) ) 
F=(F2-F1)/(2.0E0*DELZ) 

G2=U( 1 + 1, J,K+1)*W( I+ltU,K+l)+U< 1-1, J,K)*W( 1-1 , J, K) 

G1«=U( I~l, J»K+1) *W( I-1,'J,K+1 )+U( I+UJ ,K )*W( I+l, J,K) 

G= ( G2-G 1) / ( 4. OEO*-DELX*DELZ ) 

H=(2.0E0/DELX**2+2 .0E0/DELXW+2 .OEO/ ( DE L Y( J )*DEL Y( J ) ) ) 

APHI= ( 1. OEO— OMEGA)* PHI ( I , J » K ) +OMEGA* ( A+ B+C + C+2 . 0 EO* ( E+F + G ) )/H 
RESID=ABS(APHI-PHI( I,J,K)) 

ERR=MAX1 ( ERRtRES ID) 

PHI ( I , J,K)=APHI 
42 CONTINUE 



CALCULATION OF PRESSURE FROM 3D POISSON'S EQUATION 

CALCULATION ON A CORNER 

J-I AND K+l TERMS HAVE BEEN ADJUSTED 

DC 43 I =2» LEI 

J = 1 



K = NE 

Al= (PHI ( 1 + 1, J,K)+PHI( I- It J,K)) / D E L X* £ 2 
A2=(PHI (I»U»K)+PHI (I, J.K-1 ) )/DELZ**2 

A 3=2. OE 0* ( PHI (I , J+ 1,K) / DELY ( J) + PHI LI , J + l , K ) /DELY ( J ) ) / 
1 ( DELY ( J ) + DE.LY ( J ) ) 



A =A 1+A2 + A3 
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3 = (U( I + 1 » J*K) **2-2 .0E0*U( I, J,K J«2+U( 1-1 , J , K )**2 ) / DELX+ *2 
C 1= V( I , J+l t K)«^2/(JE:LY( J) 

C2 = (V ( I , J, K )-'2 )>•- ( 1 . OEO/DEL Y( J ) + 1 .Q/DELY ( J ) ) 

C 3= V( I ,J+1,K)<' Y2/DELYI J) 

C -=-C3 

C = 2 .0 EO* (C1-C2+C3 ) / ( D ELY ( J ) +DEL Y ( J 1 ) 

D = (W(I»J»Kj^? t 2-2.0E0' i W(I»J»K)<‘>2+W(I , J, K— 1 )**2 ) / D EL Z* » 2 
E^A=U( I+i, J + 1,K)*V(I + 1 ( J+1«K)/DELY(J)**2 

E2E=U(I+1»J4K)^V( I + lfJtK)* ( t. OEO/DEL Y(J >* *2- 1 . OEO / DE L Y ( J)t *2) 
E2C=U( 1 + 1, J+1,KH V(I + 1,.J+1,K)/D£LY( J)**2 

EZC=-E2C 

F2«D=LY ( JI”DELY tJJ* ( E2 A+E2B-E2C )/ (DELY ( J ) +CELY ( J ) ) 

E1A=U( I-i, J + 1,K)*V( 1-1 t J+l,K)/DcLYfJ)*r2 

Eie=U( I-I, J,K)A\n 1-1, J,K)M 1.0 EO/DELY( J ) ** 2-1 . OE O/DEL Y( J)*v.'2) 
ElC = U( 1-1 * J+l tK)*V(I-l » J+ 1 » K) / DELY (J)* *2 
E 1C=-E 1C 

E 1 =0t LY ( J ) + OELY ( J ) '' ( E 1 A + El B-El C ) / ( DELY ( J ) +DEL Y ( J ) ) 

= -( F2-E 1) /( 2. OE 0 f D L L X ) 

F 2 A=V ( !,J + 1,K)»’W( I , J+l,KJ/DcLY( J )*•*.* 2 

F2c=V(I ,J,K)'+W(I , J,K)M1 .OEO/ DELY (J)*+2-1.0E0/ DELY (J )**2) 

F2C=V( I , J + l, K)* W( 1 , J + l ,K) /DFLY( J)*s 2 

F2C=-F2C 

F2 = 0ELY(J)'-DELY(J)'- ( F 2 A+ F2 6- c 2 C ) / ( DELY ( J ) + C EL Y ( J ) ) 

F1A=V( I , J+l,K-l)*to( I , J + 1,K-1) /DELY( J)*v 2 

FI e=V( I , J ,K-1 )*W( I , J, K-l)*( 1.0c0/DELY( J ) * * 2~1 . OE O/DEL Y( J 1***2) 
F1C = V( I ,J+1,K-1)’*W(I , J+1,K-1)/DELY(J)*. 32 
F 1C=-F 1C 

Fl = DELY(J)*DtLY(J) i: (rlA+F13-FlC)/ (DELY ( J )+DELY(J) 1 
F -.( F2-F 1) / (2. 0E0RDELZ1 

G2=U( 1+1, J,K)*W( 1+1, J ,K)+U( I- It J,K-1) *W( 1-1 ,J,K-1) 

G1=U( 1-1, J ,K) +W( 1-1 , J, K)+U( 1+1 , J, K-l )* W ( I +1, J,K-1 ) 

G=(G2-G1)/(4.OEO*DELX*-0ELZ) 

H=(2.0cC/DtLX«- : f2+2 .0 EC/ DELX**2+ 2 . OEO/ ( DEL Y ( J )*DfcL Y( J ) 1 1 

APHI = < 1. OEO- OMEGA) - PH I (I , J , K) + OMEGA* (A+B + C+C+2 .0 E0-- ( E+F+G) )/H 

RES ID=AdS ( APHI-PHI ( I, J,K)) 

E RP.-MAX 1 ( E RR , RES ID) 

PHI(I,J,K)=APHI 

CONTINUE 

CALCULATION OF PRESSURE FROM 3D POISSON'S ECUATION 

CALCULATION ON A CORNER 

J+l AND K— 1 TERMS HAVE BEEN ADJUSTED 

DO 44 1=2, LEI 

J = ME 

K = 1 

Al= (PHI ( 1 + 1, J,K)+PHI( 1-1, J,K) ) /DELX4-C2 
A2=(PHI ( I t J,K+l )+PHl(I,J,K) )/DELZ**2 

A 3=2. OE 0*( PHI ( I , J ,K) /DELY ( J) +PHI (I , J-l , K)/DELY( J-l ) )/ 
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1 ( DfcLY( J ) + DEL Y( J-l ) ) 

A= A1+A2+A3 

B=(U(I+1,J,K )*$2-2 .0 EO~*U ( I , J-, K )**2+U< I - i, J , K )*•* 2 ) /DELX*'*2 
C 2 = ( V( I ,J,K)**2)*(1.0EO/DELY(J-1)+1.0/D£LY(J) ) 

C3=V( I , J-1,K)*42/DELY( J-l) 

C =2 . 0E0^ (C1-C2+C3)/ ( DELY ( J ) +DEL Y ( J-l ) ) 

D = ( W( I , J,K+1) **2-2. 0E0* W(I , J,K)»*r-2 + W( I » J , K )=**2 )/ DELZ*^2 
E2A=U(I+1, J,K)+V(I+1, J,K)/DELY< J)**2 

E2B-UU+1, J,K)*V( 1+1 , J t K)?(r.O EO/DELY (J-l )**2-l . 0 EG/ DEL Y ( J )**2) 

E2C=U( 1 + 1. J-1,K)*V( 1+1 ,J-1 ,K)/DELY( J-l)»* 2 

£2=CELY ( J )*DELY{ J-l )*( E2A+E2B-E2C )/ ( DEL Y ( J ) +DEL Y( J-l) ) 

E1A = U( 1-1, J,K)*V(I-1, J*K) /DELY ( J )**2 

E1B=U( 1-1, J,KHV( 1-1, J,K)fc( i.OEO/DELY( J- 1 )* *'*2-1 . OE O/DE LY { J ) ** 2 ) 
E1C = U( I — 1 , J-l tK)*V ( 1-1 , J-l, K)/ DELY (J-l)** 2 
E l=DEl.Y ( J ) ? : DE L Y( J-l ) + ( E 1A+ El B-E 1 C) / ( DELY ( J ) +DE LY ( J-l ) ) 

E= ( E2— El) / ( 2 .OEO*DELX ) 

F 2A = V ( I ,J,K+1)*>W(I » J » K+l ) / DtLY ( J )**2 

F2B=V( I ,J,K+1J*W( I ,J,K+1)*( i.OEO/DELY( J- 1 2-1 . 0 EO/ DEL Y ( J )4*2 ) 

F2C=V ( I . J-cl , K+l f*W( I, J-l, K+l)/ DELY (J-l) 2 

F2=DELY( J)*DELY( J-l ) ^ ( F2A+ F2B-F2 C) / ( DELY ( J ) +DELY ( J -1 )) 

F1A=V( I , J,.K)*W( I, J,K)/DELY( JM-*2 

F 1B = V ( I ,J,K) *-W( I , J ,K)* ( 1.0 EO/DELY (J-l )»*2-l .OEO/DELY( J )**2) 

F 1C = V( I » J- 1 » K ) ?• W ( I , J-l , K ) / DELY ( J-l )*'*2 

F1=DELY ( J )* DELY (J-l )* ( F 1A + F 1B-F 1C ) /( DEL Y( J ) +DE LY ( J- 1 ) ) 

F = (F2-F1) / (2. OEO* DELZ ) 

G2=U( 1+1, J,K+1)*W( 1 + 1, J,K+1)+U( I — 1,J,K)%>W(I— 1»J,K) 

G1=U( I-l» J,K + 1 )*W( 1-1, J,K+1)+U( 1 + 1, J,K)*W( I+1,J,K) 

G = ( G2— G 1 ) / ( 4. OE O'^DE LX* DELZ) 

H= ( 2 . Ot O/DEL X'^* 1 2+ 2 . OtO/Dt LX4»< 2+ 2. 0E0/< DEL Y< J)* DELY ( J-l ) ) ) 
APHl=(l .OEO-OMEGA )*PHI ( I, J , K )+0MEGA*(A+B+C + D+2.0E0HE+F + G) )/H 
RESID=A8S( A PHI-PHI (I , J,K) ) 

ERR=MAX 1( ERR, RES ID) 

PHI (I , J,K) =A PH I 
44 CONTINUE 
C 

C CALCULATION OF PRESSURE FROM 3D POISSON'S EQUATION 

C CALCULATION ON A CORNER 

C J+l AND K+l TERMS HAVE BEEN ADJUSTED 

DO 43 I =2, LEI 
J = M E 
K = NE 

A 1=(PHI ( I,+ l , J, Kl + PHI { I -1 » J ,K) ) / DELX#*2 
A2= (PH I ( I,J,K)+PHI ( I, J,K-1) )/DELZ**2 

A 3=2. OE Op ( PHI ( I , J, K)/OELY< J )+PHI ( I , J-l ,K )/ DELY( J-l ) >/ 

1 ( CEL Y( J J+DELYU-l) ) 

A = A1+A2 +A3 

B=( U( 1 + 1, J,K) +*2~2 . OE 0* LK I , J , K) **2 +U ( I -1 , J ,10**2 )/ DELX*=*2 
C1=V( I , J,K )*-*2/DELY( J ) 
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c 



c 



r 

C 



C 2= ( V ( I » J » K ) ? *2) * ( 1.0E0/DELY! J- 1) + 1. 0/D E L Y ( J) ) 

C3=V! I , J-l ,K)**2/DELY ( J-l) 

C=2.0E0*(C1-C2+C3) /(DELY( J)+D£LY ( J-l ) ) 

D = ( W ( I » J, Ki^*2-2.0E0*W(I , J, K)**2+W! I , J , K- 1 ) **2 ) /DE LZ 3 ^ 2 
E2A=Ut 1+1 , J,K)*V ( 1+1, J,K)/DELYt J )*-*2 

E 2B=U ( I +1, J,K)*V( 1 + 1, J ,K) *> (1.0E0/DELY( J-l )**2-l . 0 E0/ DcLY !J )-.»•» 2) 

E2 C=U ( 1+1, J-1,K)*V! I+1,J-1,K)/&ELY1 J-l)**2 

E2=D£LY ( J)*DELY( J-l )*• ! E2 A+ E2 3~E2 C )/ ( DELY ( J ) +DELY { J - 1 ) ) 

E 1A=U ( I-1,.J,K)*V! I- 1 , J » K ) /DE LY! J)*+2 

E1B=U!I-1, J,K)*V( 1-1, J,K)e ( 1.0E0/DELY( J-l)**2-l. OE O/DEL Y( J ) ** 2 ) 

E1C = U( I - 1 , J-l , K) + V ! I -1 , J-l , K) / D ELY (J-l)** 2 

E1=D2LY ( J )*DELY! J- 1)* ( E1A+E 1B-E 1C ) / ( DEL Y 1 J ) +DE LY ( J-l)) 

E = ( E2-E1 )/ (2.0 EO't DELX ) 

F2A = V( I , J,K)*W! I , J,K) /DELY! J)**2 

F2 6=V ( I , J, K)*W( I, J,K)* ( 1.0E0/DELY( J- 1 ) v * 2- 1 . OE O/DE LY ( J ) ‘ * 2 ) 

F 2C=V ( 1 , J—l » K ) * W ! I ,J-1,K)/DELY (J-l )*^2 

F2=DEL Y( J )*DELY! J-l)* (F2A + F26-F2C) / (DELY (JJ+DELY (J-l) ) 

FI A=V ( I , J » K-l )»*W ( I»J»K-1 )/DElY! J )**2 

F 1B=V ( I , J ,K-1 )*W( I , J,K-1 )*! 1.0E0/DELY( J-l )**2~1 .OEO/DEL Y( J H*2) 

F1C=V! I , J-1,K-1)*W! I , J-l, K-l) /DELY( J-l)** 2 

Fl = DELY (J )*DELY( J-l )* ( F1A + F 19-F 1C )/ ! DEL Y( J)+DELY{ J-l)) 

F = (F2-Fl)/(2.0E0 iI DELZ) 

G2=U!I+1,J,K)*W< I+l,JtK)+U(I-l,JtK-l)*W{I-l,J,K-l) 

G 1 = U! I — 1 , J » K) * W ! I — 1 , J , K ) +U ( 1+1, J,K-1)**{ I+1»J»K-1) 

G- ! G2-G 1) /( 4.0E0*QELX*DELZ) 

H=12.0E0/DELX**2+2 .0 EO / DEL X*? 2+ 2 . OEO/ ( DE L Y ( J ) *DEL Y { J- 1 ) ) ) 

A PHI = ( 1.0E0-0MEGA)*PHI ( I , J , K) + OMEG A* ( A+ B + C+ C+2 .0 E0-> ( E+F + G ) )/H 
R£SID=ABS( APHI-PHI ! I , J ,K) ) 

ERR=MAX1 ! ERR , RES ID) ~ 

PHI ( I , J ,K) =APHI 

45 CONTINUE 

DO 46 J = 1,ME 
DC 46 K = 1 • N E 

PHI (1, J,K) =3.0E0*PHI (2 , J,K)-3.0E0+PHI(3f J*K) + PHI(4» J»K ) 

46 CONTINUE 



PHI (LE , J,K) =3.0E04PHI ( LEI , J , K ) -3 .0 EO*P H I { L E2, J , K ) +PH I ( L E 3, J , K) 

47 CONTINUE 

IP(ERR.GT.ERRTOL) GO TC 200 

IF(IC.EQ.O) GO TO 150 
WRITE! 6,313) 

813 FORMAT! 1H0,20X, 'GUT OF CALCULATION POISSON') 

IF! ITER.EQ.l ) GO TO 49 



62 



IT ER=1 
GO TO 18 

49 T = T+DEL T 

WR ITE (6 *51 ) T 

51 FORMAT! 1H0,20X,' TI ME =',E15.5) 

DO 52 J = 1 * 50 

I=LE+1-J 

WRITE! 6,50) U(LP1,J),U(LP2,J),U!J,1),U( 

50 FORMAT ( 1H0,4E20.5) 

52 CONTINUE 
C-0 TO 14 
END 

C /* CARD MUST APPPEAR HERE 

//GO .FT06F001 DD S YSOUT = A , SP ACE = ! CYL , ( 5 , 1 ) ) 
//GO. FT16F001 DD UN IT = SY SDA, SPA CE= ( C YL , ( 5, 1) , 
// DCB = (RECFM = VS,LRECL=3504 ,BLKS IZE=3 508 ) ,DI 
// GO . FT 16F00 1 DD UN I T= SY SDA , SPACE =( C YL , ( 5 , 1) , 
// CCB= (PECFM=VS,LRECL=3504, BLKS IZE=3508), DI 
/ /GO . FT 19F 00 1 DD UNI T = SY SDA , SPACE = ( CY L , ( 5 , 1 ) , 
// CCB=(RECFM=FB»LRECL=14,BLKSIZE=3500),DISP 
//GG.FT20F001 DD UNIT=SYS DA, S PACE= ( CYL, ( 5, 1 ) , 
// DC8 = ( RECF M=FB, LRECL= 14, BLKS I ZE=35 00 ) ,DISP 



1 , 1 ) 



RLSE) , 

SP=NEW, DSN = KORE AN 1 
RLSE) , 

SP=N EW, DSN = KOREA N3 
RLSE), 

=NEW,DSN=K0REAN4 

RLSE), 

=NEW,DSN=K0REAN5 
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APPENDIX E 

I. FLOW CHART FOR 2-D PROBLEM 
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nonoonoononnnoonnooonnnooooooooo 



APPENDIX E 



II. LISTING OF CO.MPUTFR PROGRAM FOR 2-D CASE 

THE FOLLOWING FORTRAN CODED FINITE DIFFERENCE COMPUTER PROGRAM 
SOLVES THE FULL TIME DEPENDENT NAVIER STOKES EQUATIONS FOR THE 
TwG DIMENSIONAL FLOW ABOUT A FINITE, INFINITELY THIN FLAT PLATE 
IMPULSIVELY STARTED IN ITS OWN PLANE. "LEAP-FROG TIME-WISE INTEG- 
RATION IS INCLUDED AS AN OPTION. EXTENSIVE DISC WRITING IS EMPLOY- 
ED TO SAVE COMPUTER CORE SPACE. ADDITIONALLY, THE EQUATIONS 
THEMSELVES ARE APPLIED TO THE BOUNDARIES OF THE CALCULATION REGION 
IN LIEU OF EXTRAPOLATING THE VARIABLES CALCULATED WITHIN THE 
REGION TO THE BOUNDARIES. 

THE FOLLOWING PARAMETERS MUST BE SPECIFIED 
ICALC LE "ME LP1 LP2 RE DE LX DELY OMEGA ERRTGL PE QE 



t 

4 

A 

* 

% 

* 

f 

*V 

* 

* 

4 

* 

t 

* 

* 

¥ 

* 

\ 

* 



ICALC SPECIFES THE USE OF “LEAP-FROG* 1 IN T 
LE IS THE LENGTH 0^ THE REGION IN DELX STEPS 

ME IS THE HEIGHT OF THE REGION IN DELY STEPS 

LP1 IS THE START CF THE PLATE 
L P 2 1$^ THE END OF THE PLATE 
RE IS TH E REYNOLDS NUMBER BASED ON U AND L 
DELX IS THE GRID SPACING IN X 

DELY IS THE GRID SPACING IN Y 

OMEGA IS THE RELAXATION PARAMETER IN POISSONS 
EPRTCL IS THE ERROR TOLERANCE IN POISSONS 
PE IS THE FREESTREAM PRESSURE IN PSF 
QE IS THE FREESTREAM DYNAMIC PRESSURE IN PSF 
NOTE THAT DELY IS AN ARRAY DELY (ME) 

ALL VARIABLES ARE DIMENSIONED F(LE,ME) 



* * Vf * * $•*;*$** + ** ** y ** * £ * \ ± 
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no on 



c 

DIMEMS ICN U( 100, 50) , V( 100, 5 0) 100,50) ,DELY( 50) 

n 

T = 0. 0E 0 
C 

L cl =L E- 1 
L c 2=LE- 2 
L E3=L E- 3 

r 

Mcl=ME- 1 
ME2=ME-2 
ME3=ME-3 
C 

REI = 1.0E0/RF 
C 

PHI E = PE /( 2 . 0 E 0 * QE ) 

C 

A = 1 .0E0/DELX 
tJ = 1.0E0/DELY( 1) 

DELT= 1 .0E0/ < A+0. 1 E0 >E + 2 . 0E O' R E I HA *-*2+3 *>2) ) 

C 

DEL T1 =DE IT 
DELT2=DELT J/2.0E0 

C 

IC=0 

C 

C INITIALIZE THE VARIABLES WITHIN THE CALCULATION REGION 

DO 2 1=1, LE 
DO 2 J= 1, ME 
U ( I ,J) = 1 . 0 EO 
V( I , J ) = 0. 0E0 
PHI U, J ) = P HI E 
2 COM INUE 

SET PLATE TO ZERO VELOCITY 
DC 3 I = LP1,LP2 
U( IJ ) = C.OEO 
. 3 CONTINUE 

CALCJLAT= v FROM CONTINUITY EQUATION 
DC 4 J = 1 * ME1 
DC 4 I =2, LEI 
C1=U( 1 + 1, J + l)+U( 1 + 1, J ) 

C2=U( 1-1, J+l) +U( I -1 , J ) 

V(I,J+1) = V(I,J)-DELY(J) MC1-C2) /< 4.0E0> DE LX) 

4 CONTINUE 
C 

C SET UPSTREAM VALUES FOR V VELOCITY 
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on o o on o oo 



OG 5 J= 1 , ME 
V! l t J )=0.0E0 

5 CGNTINUE 

EXTRAPOLATE DOWNSTREAM VALUES FOR V VELOCITY 
DC 6 J=1,ME 

V(LE,J) = 3.0E0*V( LEI, J ) -3. 0E0” V ( L E2 , J ) + V ( LF3»J) 

6 CCNTINUE 

WRITE! 6,810) 

610 FORMAT ( 1H0, 20X, • INTO IC POISSONS') 

CALCULATE INITIAL P FIELD FROM POISSONS ECUAT ION 
GO TO 200 

150 IC=i 

WRITE (6 ,811 ) 

811 FORMAT! 1H0,20X, • OUT OF IC POISSONS') 

I T ER= 1 
14 CONTINUE 

IF! ICALC.EQ.O) GO TO 18 
I TER=0 

REWIND 16 

STORE VALUES OF PRESENT TIME ON TAPE 
WR I TE ( 16) U 

18 CONTINUE 
REWIND 19 

CALCULATE TIME DERIVATIVE OF U FROM MOMENTUM EQUATION 
DO 15 I =2, LEI 
CO 19 J = 2» ME 1 

A = ( U( 1 + 1, JM *2-U( 1-1, J >**2 )/(2 .0 EO-^DELX f 
6 1=U( I , J + 1 ) ' V ( I , J+l) /DELY! J} + ‘ 2 

B2 = U ( I, J )tv( I, J ) M l.OEO/DELY! J- 1 >**2-1. OE O/DE LY( J) **2) 
B3 = U(I , J-l )* V < I, J-l )/ DELY (J-l )■«••■* 2 

B=CELY( J jr-DELY! J- 1 ) ( B 1+B2-B3) / ( DE LY( J ) + DE LY ( J-l ) ) 

C = ( PH I ( I +1 ,J)-PHI(I — 1,J) ) / ( 2.0E0* DELX ) 
D2A=V(I+1,J+1)/DELY( J)*F2 

D2B=V ( 1+1, J )*( 1. OEO/DEL Y( J- 1) R ' 2-1.0E0/DELY(J) <l '-2) 
D2C=V(I+1, J-l)/ DELY! J-l )**2 

D2=DELY! J) FDELY! J-l)> ( D2 A+ D2B-D2 C ) / ( DELY ( J ) +DELY ( J-l ) ) 
D 1 A=V ! 1-1, J + l )/ DELY! J )*»2 

Die=V( I -1 , J) Ml .OEO/DELY! J-l )>-< 2-1 .0 EO/ DELY ( J )>* 2 ) 

D1C = V( 1-1, J-D/DELY! J-l)* '2 
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0 i = DE L Y ( J ) O E L Y ( J - 1 ) ‘ ( D1A+ D1P-D 1C ) /< DEL Y ( J) +DELY < J-l ) ) 
C=(D2-ri)/l2.0E0- DcLX) 

E 1 = U( I , J+ 1 ) / OE L Y ( J ) 

E2=U( I, J )*t 1.0EO/DELY( J-l ) + l.OE 0/DELY( J) ) 

E3=U( I, J-l l/DELY (J-l ) 

E=2.0E0* ( E1-E2+E3 ) /(DELY( J) +DFLY ( J-l ) ) 

U~EE=-( A+B+C+PE I ' ( 0-E ) ) 

WRITF(i9,9G3 ) UTEE 
90 3 FGFMAT( 1PE14.7) 

19 CONTINUE 

r 

C CALCULATE TIME DERIVATIVE OF U FROM MOMENTUM EQUATION 

CG 22 1=2, LEI 
J = 1 

A=(U( 1+ I, J )? ‘2-U( I - 1, J ) v *2) / (2. OEO'-DELX) 

B1 = U( It J+l ) A V t I , J + 1)/DELY(J K >2 

B2 = U( I , J) YV( I , Ji M 1 .OFC/DELY( J )* - 2-1 . 0 EO / C ELY ( J )* >2) 
B3=U( I , J +1 ) r - V( I,J + 1)/DELY( J ) '• »• 2 

E3=-B3 

B=DELY( J)'0£LY< J)' 1 ( Bit- B2-B3 ) / ( 0 ELY ( J ) + DELY ( J ) ) 

C= (PHI ( 1 + 1, J )-3HI ( 1-1, J ) ) /( 2.0E0- DbLX) 
D2A=V(I+1,J+1)/DELY(J)«*2 

D2B = V( 1 + 1, J ) U 1. OEO/DELY( J) . -2-1 .OEO/DcLY ( J ) ^ *2 ) 

D2C=V ( I +1, J + l ) / DEL Y ( J *2 
D2C =-D2 C 

D2= DEL Y ( J ) >DELY( J ) •» ( D 2A + D 2B-D2C ) / ( DE LY( J)+DELY(J) ) 
D1A=V( 1-1 , J+l )/DELY< J )*+2 

D 1B = V( I -1 , J) * ( 1.0E0/DELY(J)- < *2-1.0E0/DELY(J )t*2 ) 

D 1C=V ( 1-1, J+1)/DELY( J )**2 
D1C=-D1 C 

Dl=D5LY(J)DELY(J) i ( D 1A+D1 B-Dl C ) / ( DELY ( J ) +DELY ( J ) ) 

D=( D2-D1)/ ( 2 .0 p 0 *DEL X ) 

E1=U( I , J+1)/DELY(JJ 

F2=U( It J )*< l.OEO/DELYl J) + l. OEO/DELY( J) ) 

E3 =U ( I, J + l J/DELY ( J ) 

E=2.0t0'MEl-E2 + E3) / (DELY( J)+DELY (J) ) 

UTEE=-( A+B+C + PE I ; ( D-£ ) ) 

WRITE(19,903 ) UTEE 
22 CONTINUE 
C 

C CALCULATE TIME DERIVATIVE OF J FROM MOMENTUM EQUATION 

DO 23 1=2, LEI 
J = ME 

A = ( U( I + 1 , J)» *2-U( I-1,J )-*2 )/ (2 .OEO*DELX ) 

B 1=U( I , J)*V( I , J) /DELY( J) * *2 

B2=U( I, J )*V( I, J ) !■( 1.0 EO/DELY(J- 1)^12- l.OE 0/DELY(J)*T 2) 
B3=U( I , J-l) ' V( I , J-i )/DELY( J-l )«*S2 

B = DELY(J) A DELY( J-l) ( B 1+B 2-B 3) / ( DE L Y( J) + DE L Y ( J-l ) ) 
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C = C PHI ( I + 1,J)-PHI ( 1-1 , J) ) /(2.0E0' DELX) 

D2A=V ( I +1 » J ) / DEL Y ( J )'--*2 

D2B=V(I + 1 » J ) M 1 . 0E0/DELY( J- 1 ) **2-1 .0 E0/ D ELY ( J ) ‘ *2 ) 
D2C=V ( I+lt J-1)/DELY( J- 1)*?2 

C2 = DELY ( J ) : ' i D ELY ( J - 1 ) ~ ( D2A + D 28- D2C ) / ( DEL Y ( J )+DELY( J- 1) ) 
D 1 A =V ( I-1,J)/DELY( J)**2 

D1B=V{ 1-1, J ) >t 1. OF O/DEL Y( J - 1) ** 2- 1 . OE 0/D E L Y ( J ) - * 2 ) 
D1C=V( 1-1 , J-1)/DELY ( J — 1 )**Z 

0 1=DE LY ( J ) ' D E L Y ( J - 1 ) '- ( D1A + D 16-D 1C ) /( DE LY ( J ) +DELY t J — 1 ) ) 
D= ( 02— D 1 ) / (2.0E0 -DELX ) 

F 1-=U( I , J) / DEL Y( J ) 

E2=U( It J)?( l.OE O/DEL Y( J-1) + 1.0E 0/DELY( J) ) 

E3 = U ( I, J-l )/ DELY ( J-l ) 

E = 2. OE OM E1-E2 + E3) / ( DE LY ( J ) +0E LY ( J-l ) ) 

UTEE=- t A+B+C + RE Iv( D-E ) ) 

WR ITE ( 1 9 , 903 ) UTEE 
23 CONTINUE 

CC 28 J =1 » ME 
1 = 1 

UT EE= 0 . OEO 
WRITE (19,903 ) UTEE 
23 CONTINUE 

DC 29 J =1 , ME 
I = LE 

UT EE=0 .OEO 
WRITE (19, 9 03) UTEE 
29 CONTINUE 

IF ( ITER.EQ.O) DELT = DELT2 
IF< ITER .EO.l) DELT =DE L T1 
I F ( ITER .EQ.O ) GO TO 402 
IFl ICALC.EQ.O) GO TO 402 

REWIND 16 

READ VALUES OF U FROM TAPE 
R EAD( 16 ) U 

402 CONTINUE 



R EW IN D 19 

READ VALUFS OF UTEE FROM TARE AND COMPUTE NEW U 
DO 410 1=2, LEI 
DC 410 J = 2 » ME1 
READ( 19,903) UTEE 
U( I,J ) = U( I,J ) HITEE-DELT 
410 CONTINUE 
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c 



r 



r 



C 

c 

c 



DC 413 1=2 i L El 
J = 1 

PE£0(19,903) U T F E 

U( I , J) =U( I , J) +UT : E ' DELT 

413 CCNTINUE 

DC 414 1=2, LEI 
J = ME 

F EAD (19 ,903) UT h E 

U( I, J) = U( I , J ) +UTEE-DELT 

414 CONTINUE 

DC 419 J=1,M5 
1 = 1 

U ( I , J) =1 . 0 EO 

h19 ClNTINUE 



DC 415 I=LPl,LP2 
J= 1 

U ( I » J ) = 0 .0 EO 
415 CCNTINUE 

DC' 42 0 J=1,ME 
I =LE 

U(I,J) = 3.0E0 V U(LE1,J)-3.0E0‘U(LE2,J)+U(LE3,J) 
420 CCNTINUE 

DC' 30 1=2, LEI 
DO 30 J =1 , ME 1 

UEX2=( J(I + 1, J)+U( 1+1, J+l) ) /2.0E0 
UEX1=(U( I - 1, J )+U( 1-1, J + l) ) /2.0E0 
UEX = ( UEX2-UEX1 )/ (2 . 0 EO * DELX ) 

V(I,J + 1) = V(I,J)-DELY(J)*UEX 
30 CONTINUE 

DO 34 I =1 ,LF 
V( I » 1 ) = 0 .OEO 

34 CCNTINUE 

DO 35 J=1,MF 
V ( 1 , J)=O.OEO 

35 CONTINUE 



V( It J ) = 3. OEO* V( L£ 1 , JJ-3.0E0 * V( LE2, J) + V( LE3 , J) 
36 CONTINUE 
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WRITE (6,812) 

8 12 FCRMATC 1H0.20X, • INTO CALCULATION POISSON') 

C 

200 CONTINUE 
C 

E R R = 0 . 0 EO 

r 

C CALCULATION OF PRESSURE FROM 2D POISSON'S EQUATION 

CO 37 1=2, LEI 
GO 37 J =2 , ME1 

Al= (PHI ( 1 + 1, J )+PHI( 1-1, J) ) /OELXi 2 

A2 =2 . 0 E0“ (PHI ( I , J + l ) / DELY ( J ) +PH I ( I , J- 1 ) / DEL Y( J- 1 ) ) / 

1( D5tY( J)+DCLY( J- 1 ) ) 

A=A1+A2 

6=(U(I+1 , J)**2-2.0E0 "U( I, J ) - 2+U ( 1-1, J ) * > 2 ) /DELX* 2 
C 1 = V ( I , J+l)< * 2/DE L Y( J ) 

C2 = V ( I, J )* 2M 1.0E0/DELY( J- 1 ) + 1 . OEO/DE L Y( J) ) 

C 3 = V( I , J-l) *2 /DELY (J-l) 

C = 2 » 0 E 0 ' (C1-C2 + C3) /( DE L Y( J)+DELY( J-l) ) 

C2A=U( 1+1, J+l ) ■ V( 1+1, J + 1)/DELY( J )^2 

D 2B = U( 1+1 , J) * V( I <-1 , J) Ml. 050/ DELY ( J-l ) ' ' 2-1 .OEO/DELYU )-• -2 ) 
D2C=U( 1 + 1. J-1PV( 1 + 1, J-l) /DELY( J-l)>*2 

D2=DELY ( J) 'DELY ( J-l) * ( D2A + D2B-D2 C ) / ( DEL Y( J ) +DEL Y ( J-l)) 
D1A=U( 1-1, J*l) ;: VU-1* J+l) / DELY( J)**2 

C1E=J( I - 1 , J )*V( 1-1, J )* ( 1. OEO/DEL Y( J-l) ’*t 2-1. OEO/DE LY( J) -• • 2) 
D 1C=U ( I -1 , J-l )* V( 1-1 , J-l ) / DELY ( J-l )-» r ’2 

C 1= DSL Y ( J ) DEL Y( J- 1 ) ' ( D1A + D13-D 1C) / ( DELY ( J) +DELY ( J-l ) ) 

C= ( D2-D 1 ) / ( 2 .0 ED • DELX ) 

H=2. OEO/DE LX* "2+2. OEO/(DELY (J) . DELY(J-1 ) ) 

A P H 1= ( l.OEO- OMEGA > * PHI ( I , J) + OMEGA* ( A+B+C + 2. OEO- D) / H 
RESID = ABS ( APHI-PHI ( I, J ) ) 

FRR=MAX1(EPR,RESID) 

PHI ( I, J ) = APH I 
37 CONTINUE 
C 

C CALCULATION OF PRESSURE FROM 2D POISSON'S EQUATION 

DC 43 I =2 , LEI 
J = ME 

A1=(PHI ( I + 1,J )+PHI(I — 1,J) ) /DELX ^2 

A2=2. OE 0* (PHI ( I, J)/DELY( J)+PHI( I , J-l ) / DELY ( J- 1 ) )/ 

1 ( CEL Y( J ) +DEL Y ( J- 1 ) ) 

A=A1+A2 

E = ( U( 1 + 1, J) * fr 2-2. OEO U(I , J)+U(I-1 ,J)** 2 )/CELX* *2 
C1=V( I , J ) £?2/ DEL Y ( J ) 

C 2 =V( I , J) *- 2 * ( l.OEO/ DELY (J-l ) + l .OEO/DELY(J ) ) 

C 3 = V( I. J-l) s 2/DEL Y( J-l) 
C=2.0E0-(C1-C2+C3)/(DELY(J)+DELY(J-1) ) 
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DZA=U( I+1,J)&V( I + 1 , J ) /DEL Y( J ) 2 

D2 c=U ( I +1,J)*V(I+1»J)B (1.3EO/OELY(J-l)">-?2-l.JEO/DELY( J ) '■= 2) 
D2C-U( I+l.J-l) t'V(I+l,J-l) /DELY( J-l) +*2 

C2 = DEL Y (J K'DELY ( J-l ) • ( 02A+D2B-D 2C ) /( DEL Y( J)+DELY( J— 1 ) ) 

D1A = U( 1-1 , J)? V( I -1 » J ) / DELY ( J )‘ ; >2 

D1B=U( 1-1, J ) i V( !-l , J)M 1. OEO/DELY( J-l) 2-1 .OEO/ DELY ( J ) ?■ ‘2 ) 

Dl 0 = J ( 1-1, J-l ) ■ V( 1-1, J-l) /DtLY{ J-l)> >2 

01 = DELY(J)"DELY( J-l)’ ( Dl A+ D1 B-Dl C ) / ( DELY ( J ) +DELY ( J-l ) ) 

D= ( D2-D 1) / ( 2.0E0-DELX) 

H=2.0E0/DELX*< 2 + 2. OEO/ (DELY (J ) : DELY (J-l) ) 

A P H I = ( l.OEO-OMBGA) * PHI (I , J ) +CME G A* ( A+B+C +2 .OEO * D ) / H 
PcSID=ABS (APHI-PHI ( I, J ) ) 

E-R=MAX1(ERR,RESID) 

PHKI.J )=APHI 

40 CONTINUE 

r 

C C ALCJL A T I ON OF PRESSURE FRjM 2D POISSON'S EQUATION 

DC 41 I =2 , LEI 
J = 1 

A 1=( PH I ( 1 + 1, J ) +PH I ( I - 1 , J ) ) / DEL XT « 2 

A2 = 2.0E0*‘ (PHI(I,J+1)/CELY(J)+°HI(I,J+1)/0ELY(J))/ 

1( DFLY( J ) +DEL Y( J ) ) 

A= A 1+A2 

B = (U(I + l,J)M2-2.0 E0 ; U< I , J ) +U ( 1-1 , J )*-<2 ) / CELX*+2 
C 1= V ( I , J + l ) >•'- 2/DELY( J ) 

C2 = V( I, J)-»*2*<1 .0 EO/ DELY ( J ) +1.0E0/DEL Y( J ) } 

C 3= V( I , J+l)vr 2/DELY( J) 

C 3=— C3 

C=2 .0 E0-? (C1-C2+C3 ) / ( DELY ( J ) +DEL Y ( J ) ) 

D2A =U( 1 + 1, J+l)* V(I+1, J+l) /DELY (J)**2 

02 8=U ( I +1, J )*V ( I + 1, J ) M 1.0E0/DELY(J)*’ Ji 2— 1« OEO/DE LY( J) " 2) 
D2C=U(I+1, J+l ) - V( 1+1, J+l )/ DELY (J )**-2 

0 2C=-D2C 

02 = DELY (J)'fDELY(J ) ; - < D2A + D2B-D2C ) / ( DELY( J )+DELY(J) ) 

C 1 A=U ( I - 1 » J+l ) * V( 1-1 , J+l J/DELY ( J)-*‘ 2 

D1E=U( I - 1 » J ) $ V( I - 1 1 J ) *■( 1.0E0/DELY( J) M 2-1 . OEO/DEL Y ( J ) 2 ) 

Dl C-U ( 1-1, J+l ) • V( 1-1, J + l) /DELY (J )^2 
. D1C=-D1C 

D1=DELY ( J ) vDE LY( J )* (D1A+D13-0 1C ) /( DELY( J ) +DELY( J) ) 

D = (D2-D1) / <2. OEO* DELX ) 

H=2.0E0/DELX*“'2+2. OE 0 / ( DE L Y ( J) 1 DE LY( J) ) 

A PH I = ( 1 .OEO-OMEGA )'«PHI ( I , J ) +OMEGA> ( A+B + C +2 . OE 0 J D ) /H 
RESID=ABS( APHI-PHI ( I, J) ) 

ERR=MAX1( ERR. ,RESID) 

PHI ( I , J ) = APH I 

41 CONTINUE 
C 

DO 46 J = 1 , ME 
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1 = 1 

PHI (I , J) =3.0EO PHI (2 , J )-3 .0 EO > PHI (3, J ) +PHI (4, J ) 

46 CONTINUE 
C 

DC 47 J=1,ME 
I = LE 

PHI (ItJ )=3.0E0*PHI(LE1»J ) - 3 . 0 5 0 • PHI(LE2»J) + PHI ( L E 3 » J ) 

47 CONTINUE 
C 

I P(EPR.GT . ERRTGL ) GO T'J 200 
I F ( IC . EG. 0) GO TO 130 
C 

WPITE(6 ,813 ) 

613 FORMAT! 1H0,20X,’ OUT OF CALCULATION POISSON*) 
r 

IF< ITER .EQ.l ) Gu TO 49 
I TER = 1 
r 

GO TO 13 
C 

49 T = T+DEL T 
WRITE(6,51) T 

51 FORMAT! 1H0,20X, * TI ME =*,E15.5) 

r- 

DC 52 J =1 » 50 
I =LE+ 1- J 

WRITE (6, 50) U(LP1, J),U(LP2, J) , U(J,1),U( I ,1) 

50 FORMAT ( 1H0 »4E20. 5 ) 

52 CO'NT INUE 
C 

GC TO 14 
END 

C /■• CARD MUST APPPEAR HERE 

//G3.FTG6F001 DD S YSOUT = A , S PACE = ( CY L , ( 5 , 1 ) ) 

//GO .FT16F00 1 DD UN IT = SY SDA , SPACE = ( C YL , ( 5 1 1) t R L SE ) , 

// CCB=(RECPM=VS, LRECL = 3504t 6LKS IZE= 3 508 ) , D I S P = N EW , DSN = KOP E A N 1 
//GO .FT19F001 DD UNI T = SY SD A, SPACE = ( CYL , ( 5 , 1 ) , R LS E ) , 

// OCB= (RECFM=FB, LRcCL= 14, BLKS I Z E= 3500) , D I SP =NE W ,D SN=K0REAN4 



73 



INITIAL DISTRIBUTION LIST 



No. 



1. Defense Documentation Center 
Cameron Station 
Alexandria, Virginia 22314 

2 . Library- 

Naval Postgraduate School 
Monterey, California 93940 

3. Dr. John M. Wozencraft 
Dean of Research 

Naval Postgraduate School 
Monterey, California 93940 

4. Dr. Richard W. Bell 

Chairman, Department of Aeronautics 
Naval Postgraduate School 
Monterey, California 93940 

5. Professor Allen E. Fuhs 
Department of Aeronautics 
Naval Postgraduate School 
Monterey, California 93940 

6. Professor Gustave Hoke ns on 

Assistant Professor, Department of Aeronautics 
Naval Postgraduate School 
Monterey, California 93940 



74 



Copies 



20 



2 



2 



1 



1 



5 



UNCLASSIFIED 

Smmtv C 1 ,issi fn.it ton 



DOCUMENT CONTROL DATA • R & D 

t /.tssy/n-.-ifion of title, tuul\ »>/ abstract stud index in,* .tnnot..tinn must ttc unt'-ri'd iv/irn th» ‘ oivr.,// r»p.,rt i . . /., 

MjiN & ting activity ( I’nrpor.ifr /iiifhor ^ 

Naval Postgraduate School 
Monterey, California 939^0 



RU'ORI SCCUKITY CU'.SIll- /. : IO’. 

UNCLASSIFIED 



2I>. GHOUP 



None 



REPORT TITLE 



Numerical Experiments with the Two- and Three-Dimensional Unsteady Navier-Stokes 
Equat ions 

Ol script i vi NOTES (*Tyg i«* <>/ r«*|»i»r/ inclusive tlntcs) 

Technical Report, 1973 

authORiSI (First name, middte initial, last name) 

Gustave J. Hokenson 



REPOR T DATE 



April 1973 



7a. TOTAL NO. OF PAGES 



81 



7 b. NO OF REPS 

11 



a. CONTRACT OR GRANT NO. 



!>. PROJECT NO. 



9a. ORIGINATOR'S REPORT NUMBER(S) 



NPS-57Hw7304lA 



9b. OTHER REPORT NO(S) (Any other numbers that mav be as stoned 
this report) 



). DISTRIBUTION STATEMENT 



This document has been approved for public release and sale; its distribution 
is unlimited. 




ABSTRACT 



The two- and three-dimensional unsteady Navier-Stokes equations are solved 
numerically for the flow field about an impulsively started flat plate. In 
attempting to obtain an exact time dependent solution, several significant 
results were observed. First, with regard to the formulation of the differential 
equations themselves, it appears that Poisson T s equation for the pressure field 
is a fundamental equation in as much as it allows us to solve for pressure moat 
exactly at any given time. Secondly, the difference equations must be carefully 
and consistently formulated. In this research, a non-uniform lateral grid, a 
unique interpretation of the continuity equation, and "leap frog" integration in 
time proved to be valuable techniques in obtaining an exact solution. 



D 



FORM 



T NOV 65 

N 01 01 - 807-681 1 



1473 



(PAGE 1) 



UNCLASSIFIED 



75 



Security Classification 



A-31408 



UNCLASSIFIED 



S< . untv C'l.issilu'ution 



key wo nos 



Numerical Analysis 
Fluid Flow 
Navier -Stokes Equation 



ROLE 



DD , F , r „,1473 ( BACK ) 



S/N 0101 “807-6821 



76 



UNCLASSIFIED 



Security Classification 



■ .. 1 4 0 :• 



DUDLEY KNOX LIBRARY 




3 2768 00391411 0 



